{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "329e7983-8728-48dc-b21a-0efadd3c4f3f",
   "metadata": {},
   "source": [
    "# 4. Predicting the whole Kano and Yoshii River\n",
    "*Written by Men Vuthy, 2022*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39d7aa33",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3c60f3c",
   "metadata": {},
   "source": [
    "**Import modules**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1b59693e-f39b-4635-956c-30db8ab0f140",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "np.random.seed(0)\n",
    "\n",
    "import rasterio\n",
    "import geopandas as gpd\n",
    "\n",
    "# Import scikit-learn modules\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, classification_report\n",
    "import joblib\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3396a996-756c-454f-a3f0-94edb85c651d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input classified data of each river\n",
    "Kano_classified = pd.read_csv('data/kano_river/out_img/classified/kano_classified_index.csv')\n",
    "Yoshii_classified = pd.read_csv('data/yoshii_river/out_img/classified/yoshii_classified_index.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9498c8d6-45f5-43e5-aa73-eb7567a3db50",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['B1', 'G1', 'R1', 'NIR1', 'NDVI1', 'NDWI1', 'BSI1', 'B2', 'G2', 'R2',\n",
       "       'NIR2', 'NDVI2', 'NDWI2', 'BSI2', 'B3', 'G3', 'R3', 'NIR3', 'NDVI3',\n",
       "       'NDWI3', 'BSI3', 'B4', 'G4', 'R4', 'NIR4', 'NDVI4', 'NDWI4', 'BSI4'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get index of dataframe\n",
    "kano_index = Kano_classified.iloc[:,0]\n",
    "yoshii_index = Yoshii_classified.iloc[:,0]\n",
    "\n",
    "# Create a list of the feature column's names\n",
    "features = Kano_classified.columns[1:29]\n",
    "\n",
    "features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "34ddabdc-398a-4366-ac18-3adcade4df5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load, no need to initialize the loaded_rf\n",
    "rfc_model = joblib.load(\"./random_forest.joblib\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6de62334-4254-44f1-97db-edd73b4a3dfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the trained Classifier to the prediction dataframe\n",
    "preds_kano = rfc_model.predict(Kano_classified[features])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a55e1c2f-1056-4578-b8db-5c9e47d509a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our classification accuracy is: 96.09589654544772%\n"
     ]
    }
   ],
   "source": [
    "# View accuracy classification (cross-validation) score\n",
    "print('Our classification accuracy is: {cv}%'.format(cv=accuracy_score(Kano_classified['label'], preds_kano)* 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3a4c9c-29ff-468e-b4dc-3fe2a46ef5fa",
   "metadata": {},
   "source": [
    "*Visualizing confusion matrix*data/kano_river/out_img/predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7ca58219-944a-4131-8348-94d73b67fa81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get and reshape confusion matrix data\n",
    "Matrix = confusion_matrix(Kano_classified['label'], preds_kano)\n",
    "matrix = Matrix.astype('float') / Matrix.sum(axis=1)[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a5d51dcd-5232-4486-a89f-194b2fae7264",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxEAAAFZCAYAAADn86JOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABYsUlEQVR4nO3dfXwU930v+s8PuQ1I7INx0/RYIxKL3JiVcHyNwZJo2tPWEo7be9qeIOT7OmlizKPTNAHCo21AIGSDhGSE08QStmwn7e0pEsppb0+Og1e0PW2KJPOQppa0OA1yglZu01zDPgB2TsP53T9mdjS7mn3U7s6s9vPmtS+085uH7/xmdna+8/vNrJBSgoiIiIiIKFXzrA6AiIiIiIgKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIpsSQiwXQniFEBe1/3fPcj7eNKapFEJcyWR5ceZ3RQjRbzK8TQgxq+dMZxJrJnViMo/r2jz6tf97Mp2XybyzVveGeV437EuRlzvby0khjqztC4nqSQjhFkJczCRGbfrrxvrR9pkruayzRPtUtj+TRESzdYfVARDRTEKI5QD6ATRIKSe0YZUZzq4NwJbIfFKhjbskw+XFYxZ/PYBAsgmFEJXx4s8w1rTrJI61UsoAAAghNgsheqSUW2Y5z1x6OBJvLiTaTjEy3hesoH32Ip/HQI4XZ7pP5egzSUSUMbZEENnTS4g5yZ3FCW+qJ3a5NiiEaIy80RKlCylOO+PK9Szlok4GAazI8jwLTarbaTb7Ql5pLQ9eqCf3+f4ccZ8iIttiEkFkT5VSykGzAu3qZKQ7Sr02zK11VenR/m/Thu8GUKmNuzy2i0eke4TWVeKiobtL7HgpLzOBUwAeM7x/DDEnnYZlXIx0G9HmG+l+tFtbbuTvfmOsWiybtb+Xm3Vnia0Ts/WLXUaS9YqcaLYBMHY/MVuXuHVmGNYDYJFheKK6bzP8vzuyvGTxxsQeb/5R66/9bdw/YveZqO2UZLGp7Asz4kpST1HxpVMHSZwF0COlvBQTX7rb13R94ondp2L28/7Ivqu9j3yOY7dRWvsxEVFapJR88cWXjV5Qu3pciVO2HEC/4f1FAG7tJQG4teHXDeNcMfztBnAxtgzqyUqjYRx9vEyWaRL3FZNYrsSbDsBmqCdu+jJj1kEC2B1nnS5qMV+MxBYvngTr9zHjMuLM4zrUK9RebdzlccbT1yVenUHtyuM1zjvNupeG7ec1iyUmXi+A3SnMf7chDq9h/+yJ3Wdit9Ns9oUEcSWqp9j43KnEk8L27U80nxS3r+n6pLNPIfoz2QigzTjvBHWQcD/miy+++Mr0xXsiiGxGSjkhhFgUp3gL1Ku4EYNQT6wGAUzI6f7a19JcrBdAjxBiJdSTD+P02VzmaaF2Y5nQpo8i1L7n9QAakswnIKVsj1O2FupJ2iaZWv91s/X7jSTL0JclpQxoV3krAehXqxOsi1mdNSD6SnxkeKp1b6zPCRiu0JvFa4ixJ8H8jev/GIBFhqvZlQD2IHqfCSA9ifaFeOu9Eub1ZBZfNkxIKbdoV/43SylPRgrS3L7x1ue0yTLj7lMRUsrTWivHHqjrHmnZMauDVPZjIqK0sTsTkT1NpNLlQRM5YQxkujCpdp2KnAyl0h0m02X2QD2hMuu+ErmiegGGbkFxJEpY3Fr5bE4kI/NI1R6oV+YBJF2XwCziAuLUfYoJUzrzj13/I1LKtdrrwQz2mVhx94UkccUTFV+iEbXuRpGuRYn2kz3a/5sAtBm6LWVj+yZbn6h9ysSEFnujlDKSjJjVQboXFIiIUsIkgsieNkG9yuuODDCcuBif/hO5apwOtzY/N7QTbaE9VUdKuQfqiZHbMH42lglAvzl8EYB6OfOejxUATkm17/mME7uYfu6BBIt5CcCDAFYa+40nYLZ+f5tkGVG09RoU2v0YSLIuJrxQW1Ai6xk5wcxa3ceRaP4Bw/BTxvG0+yGi9hnjfQGpLDjJvhAvrnj1NCO+JMveIqVs0F5Jb5bWErQjUPctIP3tm/Z2NNmnYvVATTQiLRXx6iCQQnxERGljdyYiG5JSXhJCrAXQb+jadEpK2S6ibxg+onV/cqc434AQ4pJQfx/hEqZPQOq15QFq95KAYZrB2SzTRKSvdqw+ABeFEA3a8o1XUAcBnBVCXMD01eEZtK4ce7T13KRNk/CxpmbrB+BHqa+Org3qSe7JJOsSL4a1WgwXoG6DXNS92XKTzl/bH/vF9O9qeAEEjPuMVuf6dpKpPerWdF+IFxfUq+9m9WQW38nY+c6G9tm7oiWmmWxfs/VJxrhPxRqEmpys1ZaR8zogIjISUs7qd56IiIiIiKjIsDsTERERERGlhUkEERERERGlhUkEZURM/7CVV+snfMXwPt6NgLHzqIz8SNJsxslELuaZ4nKvx94srdWdO/5UGS9rs/bDUzN++EvrOz3jyS/a8MbY4SbjRf0YnWG4vr2SbTutLryRftza40ZnzMcuxPQPiC2PGZ7S/p7Gciz7XKQrmzHkq35j5p2XetQ+48l+jLHo2WW/NmPF/klkd0wiKCNSyvbI002gPuu8x/C0k5Ru5tOe7LJktuPkSrInvGRp/v0AGrL4eE6jBu0xj4+ZlG0CsNm4jkJ9pKzb8LjItBm3V4rbLvI4ygaoN6r2pDFtxgwJVuR1PZI8CSGkISGO/FrzcqhPAloL4CnDfHZDvck23nLS3ofs/rlIlZj+Ne2LIvoXpy2tX+M0+ahHbXkTUH8gjhJId3uI6F8BN/6Cd6V2ccL0Yok2TqOY/sXx+kTzzHT/JJrr+HQmovj6oT4qNOu0lgcv1JPoVJ7Sksn83fFaOLQn6RyB+nScyLP+jX9bYRDRj8HMGS3RPQnoddVvSJ4uaUmNUSXUJxBFns0fmS6QJAHM2T5kZ9pJ2RUp5R6tvi4CuFMrtrp+871NtkD7hW8hxHLtsbA0S9p+slZK2aDtK2cxvV37AcR9Kps2/lOR39LQWj+WJJhnpvsn0ZzGlgjKOq2ri1eoXZ4iVxojV3YuCsPz5CNdYiJ/G5qM21IdRytrMyyjx6yJ2TBdDww/9BQntjYAyyPrEW+8WTgLtfUm6oQiUT3FWW/jVTP9apr2xVYJ9bn2ewzjuw3jtEP9hdtGbR1Pa1+SM+apXdmLXFWOPELSncr2SqUytLjaoP1wV8x8+mOuMka6S+2OxKONP2O/S1FUHcVxCer+YHy+f9QvGJusU9Q+ZOPPRbr7nOnnKMZEpG4Mj2ZN1Gow6/rVhsXuE1H7bZxtkpV6TCDyY3CnEdMqKKJbxHabDYv9HBn2/5T2pzjz7BfRV99ndCEy+8wniDf2eGEWW9S2MauomO1hdsyJqleoF2Iix7tr2jT1UB+5G2kJM9tWTYj+nY4JoR5jTOeJDPZPoqIgpeSLr1m9oJ787Ta8dwOQxmGGss1QT54j412Mmcatvb+exjiVUK8kR2KpN1luPQCv4f31RLFp7y/GWd+o8dKsq+tQv6T6480/Tj2ZrffyyHpH4o2MY1iWcbssj41bG3Yd6lXjuPPU6rXREE+q2yvROkbqwqvNa3nMPhSZTyOANmN82v9ew/aP/N6A6X6XIAZ9PjFx9WvrXm8Y3qgNd2vTRepjeYL5X4xZJ1t9LjLY55J+juLM+0oe6tdsn4jab+Nsk6zWY0x8lYaYlhvrS3sfG4vZMHfMsCsx8Sbbn8zm2WhYp6jPvGG82M98vHjNjhdRsZltmzj1ZdweM7adyTq2xSy7HsBuRO9vV2Kn18bZbHjfo9WJ6Twz3T/54muuv9gSQbkSkOqVbgD6VaXNULvLJLp6GdD+jvfjTfHGcWv/v2f426gB6hdAhD5tirGlPF4KJqSUa6Fe/Yq6SpZgGWbrvQXqr9RGDEL9EoV25e4koq98PgXtSn+EVFtCLkD9wk40Ty/UK3tthrhS2V7JRO6HOI04v/wr1au4kf7kj2nxPQa1FaVfi32FVh6136Ug0tXE6KRU79N4EIZ9Rkp5WhseAPCYlPK0tvwtIvWbZu32uUh3n4v7OYpHq6O1hkG5ql+zfcJsv00k43qMYwumr2xfgtp6F2lVi+zL0MoDcYYlksr+NGOe2mfKNA6D2Loziy3uMSgmtnif10SSbbtBAI1a60W9YX0A9XgScSnF5SWcZ5Y+/0RzCpMIyhXjSXrk6vEFzDxhMwqkMN8Z40i1u8Qi7cR5pUzjxuBUY0tlPK2bQ6RZP1HXjUjXmU1QvyTdKSwjkHRlVIu0k4hLUso9AC5o3QgaAVyT5v2xA0h8MrhISjmI6fslIl0rUo0pFXswnciYiXSHaTRs3yPal3rkhBRIP5mJ7dYArd4iromZT2OpB3BK226VUv1l5nqkxlafiyztc3EJtcvTEeN+l+P6jdon4uy3iQRiB8zm+AI1+V2pdR/qh3qDtdmDDjKVyf4UMajVdb3ZOmVQdxGRE/7Yz6LZ5zWuZMvXtsseqN0RI12O9F8xTyIAwHgDt37PQ5x56ma5fxLNKUwiKFcChr9XADilnUhk/YlH2snlKak+GWptnNG80K6Gagf/yBddwtgMfXeTroOUcoucfkJV0i8y7YrWEahfWCktI0Y/om9ENvbXdUdi0sZp0/7OaJ5CiEqpPjllD9QvVXcK80qZVl+DcfovA+pJ0R6oX+qAevVTj9OQtAUyWLY+jRCiXkTfOLnIJPFq0Ial1BoV0//bGJ8dPhfpxhDvc2S2/B7E3PeT4/qdsU/E7reGhN0dO58E80+lHuNNFzCcOK+FWneRfTzSoha1HJNhgPZ5i5y4GhYTMPwdb1vGm2ePFo/pscrkM/+6yXwSHYOMscX7vMYVb9sZRVoHtPgmDK0sxpP65VpCYtQXM05lZD80m2fMtGntn0RzGZ/ORPnQB/XxnQ1IftU7E5Gb6CJfUqdjrnZCSjkohFir3bRnvFqVKLZBAGeFEBegnsBmfR2klO1CfYb88iSxmE07KNTfmYhcpTuinYyf1K58RoafhtqNIvI+7hOh4s1TqDdPRk6gJpDdVoiINqgnqWY3Kw5CPWFZq8V5SVvHyA2XptMJ9YbRB826hYjpR28aXQDwkuEE+eGYaTZDu8qr1Ytbm0+8bRW7DxlZ/rlIN4YEn6MoWj01AVghhIgM3oQc1q+UcovJPhEw7rdSfSpZom1iJm49Jtq/YNJVTttvrwntKU1a62XkszYo1adZmQ27pK3XJUwn0rFMt6XZcgDs0YbXJ6iD+pi6+9s4sZkdL9wm6x27bZJ1O4xdfiB2BMP8Il1EI44YYjImL/r2EkJExrkGQ3e7BPPMZP8kmtOElNLqGIhmRaj9Ua9IKU9qX14vQb0il/HvHdDcIIToSbEFJtX51RuvamrJ31OI6bJjB4X4ubBj/Saqx2zvX/kmhPDKmY/bpTjsuH8SWYktETQXuDF91S0ghJgAm5pJldUbHmO7RWgnDil3cckzNwrsc2HT+nUjfj0W7A21WkIUsDiMgmLT/ZPIMmyJoIJnuDro1gZNFPLVQaJs4OciO+ZiPQr1Nx4aoHVrsjoeIipMTCKIiIiIiCgtfDoTERERERGlhUkEERERERGlxU43VrNfFRERERHlmkg+ivVEg5LRubH0+vOyfrZJIsSji60OoWjJ16/i/ds3rQ6jaC0oKcMHt29ZHUZRml9Syrq3EOvfOvNLSgGAx36LLCgpY91baEFJmdUhzAm2SSKIiIiIiEgj7N1gwiSCiIiIiMhubH7nMpMIIiIiIiK7YUsEERERERGlxd45BJMIIiIiIiLbYUsEERERERGlhfdEEBERERFRWtgSQUREREREabF3DsEkgoiIiIjIdubZO4tgEkFEREREZDf2ziGYRBARERER2Q7viSAiIiIiorTYO4ew+8OjiIiIiIjIbnLaEiGEcAN4SUq5NpfLyYfuLx3BIocbE/96FXtfOWJaDgCBmyHsfeUINj36X7D2U78DAFjkcOPU3/0Vjp3uzmvMhehwcyuCwRCUinJs27E1pfJ407za+xrG3hpHR1e7Ph4AOJwO03kT0NLcilAwiPIKBdtN6sis3GzYzm27EA6HUa4oOHBonz4eoNa/2byLVbbqPNXtAEQ+G2Po6DqWhzUsDLna90/3DcB7xgsACAZDeOTR1Xhiw7q8rVehyMaxX63rQQBAKBjEaq2ueexPLlvfvTu37QYAuFxO7NeOOTu37UY4HIailOvDKEU2v7E6Zy0RWgKxGcDyXC0jXx7+Pz+FK//yYzQ99wW8F7qOB5Ysiypf86nfxsUfvoUnv/oUAOCBJcvw0ut/htXPfBarn/ksLvzzPzGBSMHwuREoixV0dLXD5XbBN+5LWh5vmlAoDP9Vvz6t98wgPNUe/QAWO28Chs8No2Kxgo6uY3Cb1v/McrNhp/sGcKDlAHp6u/XpfOM+rN+4DgcO7cPUpJ/1r8lWnae6HQD1szFp+GxQbvf9xqY16OntRk9vN6qqq5hAmMjWsV+t6xfR0/uiXtc89ieXrfo/3TeA2lU16Ohqh6faA++ZQe0zsR89vS/q86I0iAxfeZKzJEJKGZBStgMI5GoZ+dKw/Ndw6YejAICJf72K+gc+FVW+8hP348IP/gkAcP4H38eKT3xSL4skGJTc8NAwPB4PAEBRFAwPjSQtjzfNic4TeGLjOn3asdExVC+rAgBUL6vG2Oh4rlen4AwNjcDjWQoAUJTyGfVvVm42rLFpDZxOBwCgYrGCcDgMT5UHSoWiJnf+KXiqPHlcM/vKVp2nuh0AoKvzBNYbPhuU230/wnvGi6pq7vdmsnnsB6YvGgE89qciW/VfvawKvjEfQqEwfGM+eKqWRn0mlJjPBKVAiMxeecJ7IlLgLnPiWjgAAAjcCOEu551R5cbEYeUn7oe7zKmXPdX0Rbz0+p/lLdZCFg6F4XKrdedwOBAMBJOWmw3zjfvgcDqgVCj6tMYvj7HRMR7ITKh16QKg1mXAtP6jy5NNM3xuGDV1tQDUK4C7tu9CTV1NrlelYGSrzlPdDmafDcr9vg8AvS+9isamNblcjYKVrWN/xCsvvaLXNY/9yWWr/iMXh3Zt3216nBk+N8Ljf7rmicxe+Qovb0sqYIGbISxyuAEA7oVOvBe6HlU+8N3/gbWf+h10f+kIKn9lMQI3QwCAe35lMSb+9Wq+wy1YDqcDwYBad+Hw9Bd0onKzYV2dL8A/OYXDza247LuM030DaHikHt4zgzjc3Ar/5BQcDkd+V64AqHWpfnmEw2G4Tes/ujzRNC3Nrdh/cJ9+FcpT5ZnRtabYZavOU90OXZ0nMDXpR0tzK3zaZ4Nyv+/7J/1QlPJ8rEpBytaxH1DrulyZPnnlsT+5bNX/6b4BeKo96Ol9ES63S78/BVDvn9h/8Bn9M0EpKtbuTKkQQmwWQlwQQlzA5A0rQ0nIe+nvsfzj6n0QKz9xv961yWj1M5/V74no+7v/DgBY/vFlOP+D7+cv0AJXW1cLn0/tizk2OqY3lSYqNxvW0/siOrrasXXHViz1LNWvSPX0vqj3i1396dX5Wq2CUVdXA5/vMgBgdHRc76qRqDzeNMe1LjORK1HGPrZutwt+/1TO16cQZKvOU90OPb3d6Og6hm07tsJj+GwUu1zu+4C6/1ffV52PVSlI2Tr2A4Bv/DKq76uKmp7H/sSyVf+xrTzBoJpkd2ndi9kCmoFi7s4khGgDUCmEaNNutI4ipTwppVwhpVyBioW5DGVWzv7jd3GX8079CUxn//G7AIC+p9UbhVxlTnR/6Qi6v3QEp/7urxDUWiKMrRKUXO2qGgQDQf1JGrWr1GbPyNMezMrjTRMrFArjcHMrDje34pFHV/NqiInaVbUIBIL6U5RqV6ldMXZu2xW33GzYq72vYWRoBIcPtmLLhifhPeOFw+FAS3Mrdm7bhbG3xnjyqslWnae6HchcLvd9APD7eQU8kWwe+/1+f1Rd89ifXLbq/4kN6zB8bgQ7t+3G2FvjaGxaY/hMPIstG74Q1TpBKbB5S4SQUuZvaQmIRxfbI5AiJF+/ivdv37Q6jKK1oKQMH9y+ZXUYRWl+SSnr3kKsf+vMLykFAB77LbKgpIx1b6EFJWX2fnaqRjx+b0bnxvIbb+dl/fiL1UREREREdmPzVIdJBBERERGR3eTx/oZMMIkgIiIiIrIbmz9DlUkEEREREZHdsCWCiIiIiIjSYu8cgkkEEREREZHtsCWCiIiIiIjSwnsiiIiIiIgoLTZvibB5jkNERERERHbDlggiIiIiIruxd0MEkwgiIiIiItuZZ+8sgkkEEREREZHd2PyeCCYRRERERER2Y+8cgkkEEREREZHdCLZEEBERERFROphEpEi+ftXqEIragpIyq0MoavNLSq0OoWix7q3F+rcWj/3WYd1TMjbPIeyTRBARERERkWqezbMI2yQR1372b1aHULQWfeiXIZ5YanUYRUu+ehnv375pdRhFaUFJGT64fcvqMIrW/JJS1r9FIi1APPZYY0FJGeveQoXSCsTuTERERERElBYmEURERERElBYmEURERERElBab5xBMIoiIiIiI7IYtEURERERElBYmEURERERElBYBJhFERERERJQGtkQQEREREZEtCCEaAQQAVEopT6ZbHjEvVwESEREREVFmhMjslXieohEApJSD2vv6mPJ6ABNa+YQQYnm8eTGJICIiIiKymXlCZPRKYiWACe3vCQCxScIFAP1a8lAppbwUN75MV4yIiIiIiHJDCJHpa7MQ4oLhtdkwW3fMYu4yvpFSBgD0AOgH8GCi+JhEEBERERHZTKZJhJTypJRyheFlvK8hAGBRgmXWAxiUUi4BEIh0fzLDJIKIiIiIyGZycU8EgPOYbo2oBOCNKV9u6MJ0BAkSjpw+nUkI0Q810Akp5ZZcLisXjrYcQygYQrlyN764/QtJy//i9P+Lv37jbwAAoWAI9Y/8Fv5g/WdxtOUYAMDhWGg6H0pd9+OHsKjMhYmfTmJvf6dpOQAEboVMyymxw82tCAZDUCrKsW3H1pTKzYZ5zwxi+NwwAGDrjq1wOh0AgFd7X8PYW+Po6GrP0xoVlpbmVoSCQZRXKNhuUv9m5WbDdm7bhXA4jHJFwYFD+wAA3jNeDJ0bAQBsM2wTUuWy7luaWwEADqfDdN6UnWOPb9yHXdv3YKlnKQDgQMt+OJ2OuMcjmpatY/9hw75unA+P/ZnJxSNepZSnhRC7tRYHt+EGa6+UsgHASa370wSsejqTFsAmLaAZd3/b3ZvDF6Aod+O5zsNwuZx4e/ztpOW/3/i7eOHkcbxw8jiWVi/FH6z/LN4efxufW/9Z7D2wC1P+d2fMh1L3cFUdrvzbVTR9fRveuxHAAx+tiipfs+IRXPzRGJ78RjMAzCinxIbPjUBZrKCjqx0utwu+cV/ScrNh/kk/zrz+BvYf2of9h/bpX9ihUBj+q34rVq0gDJ8bRsViBR1dx+A2rf+Z5WbDTvcN4EDLAfT0duvTRbbJgUP7cMCwTUiVy7r3nvGiqtqjJxSx86bsHXsAoH71w+joakdHVzucTkfc4xFNy1b9e88MwlPtwf6YfZ3H/sxl2p0pGSllu5RyUErZbhjWoP0f0LpDDSZKIIAcJhFaAAHt7RXMvJHD1s4Pnccnqu4FANyt3I3zwxdSLv/rN/4GS7Wye6vuRblyN8KhMN71v4t7teGUvobqVbj043EAwMRP/aivqosqX3nPMlz40SgA4Pw7o1jxsWV5j7GQDQ8Nw+PxAAAURcHw0EjScrNh3jcG4XI5cbi5FVs2TLe8neg8gSc2rsvPyhSgoaEReLQrqIpSPqP+zcrNhjU2rdFPlCoWKwiHw/C+MQiny4WW5lZs2fBkHteqMOSy7kdHx1G9TL2gsWxZFcZGx/O1WgUjW8ceAPBPTqGr84R+AhvveETTslX/Y6Nj+r5evaxa39d57M9crpKIbMnXPRENAAbztKysCIdvwKV9GSx0OhAMhlIu/2bvn+L3G39Xf//2+Nt4ZucBrKxdkYfI5y53qRPXbgYBqN2V7lrojio3Jg4r71kGdymvOKUjHArD5XYCABwOB4KBYNJys2HBQBDKYgX7D+1D7aoaDJ8bgW/cB4fTAaVCye9KFRC1Ll0A1LoMmNZ/dHmyaYbPDaOmrhaBQBAVi9XuNbWravWuHaTKZd0bE4fR0XGEw+Fcr07Bydaxx+FQj/nrN61HV+cL8E/6TY9HFC1b9W9MHMZGxxAOh3nsn6WiTyKEED0AthhaJYxl+iOovvHyN3MdSlocjoUIhtSD/Y1QGC6XM6XyKf+7uFu5O2rce6vuxQsnjwNQu0FRZgK3QlhUpn5pu0udeO9GIKp84MIZrF35CLofP4TKD1cgcItf1ulwOB0IBtRkOByePkFKVG42zOV2RV2h8vl82hf6FA43t+Ky7zJO9w3kcc0Kg1qX6pd3OByG27T+o8sTTdPS3Ir9B9XuG263K+qquc93OR+rVDByWfcNjzTAe8aLluZWTE369RNdmpatY49SoejdmGpX1WB4aMT0eETRslX/DY/Uw3tmEIebW+GfnILD4eCxf5ZydGN11uQ0iRBCtAFok1JOmJUbH0H1+MbP5zKUtK2sW4kfaPcv+MYu612XkpW/Pf42qqqX6uMZ74FwuZx41/9urkOfs7xj57Bcu89h5T3L9K5NRqs7Nuj3RPS9+Xpe4yt0tXW1+hfs2OiY/sWbqNxsmMfj0Yf5/X4oioKe3hfR0dWOrTu2YqlnKRqb1uRxzQpDXV2NfnI/Ojqun/QnKo83zfHOE1i/cZ1+9c/jWaqP5/dPQVHK87JOhSKXdQ8APb3d+j0Rqz+9OufrU2iydewx9uX3X1WPPWbHI4qWrfoHgJ7eF/V7IlZ/ejWP/bNUtC0RQojdAOoB9AghvImeM2tHD9WuQDAY0p+s9JDWFenpHfsTlr/rfxcLDTduLXQ6cLTlGJ7esR/jY5ejujlRes6OD+GuhW79CUxnx4cAAH1/2AUAcC1woPvxQ+h+/BBOvfk6gu+zJSIdtatqEAwE9adr1K6qAQDs3LY7bnm8Yf6rfuzcthv+q340PFJQz1SwTO0qtdtRi16XtQDUp/3EKzcb9mrvaxgZGsHhg+r9D94zXtSuqsXkVT92btuFyat+NDzSYMEa2lcu6z4UCqOluRUtza145NHVvLHXRLaOPQ6HAzu37cbh5lY4tNYIHo+Sy1b9h0JhHG5uxWHu61lj9yRCSCnztrBErv3s3+wRSBFa9KFfhnhiafIRKSfkq5fx/u2bVodRlBaUlOGD27esDqNozS8pZf1bZH5JKQDw2GORBSVlrHsLLSgpy2Onn8xVPPebGZ0bTz79N3lZv5z+TgQREREREaUvn/c3ZIJJBBERERGRzeSza1ImmEQQEREREdmMgL2TiHz9TgQREREREc0RbIkgIiIiIrIZdmciIiIiIqK0MIkgIiIiIqK02DyHYBJBRERERGQ3bIkgIiIiIqK0MIkgIiIiIqK0MIkgIiIiIqK02DyHYBJBRERERGQ3bIkgIiIiIqK0MIkgIiIiIqK0MIkgIiIiIqK02DyHsE8SsehDv2x1CEVNvnrZ6hCK2oKSMqtDKFrzS0qtDqGosf6txWOPdVj3lAxbIlJ0/Wc/tTqEonXnhz6MWz+/YXUYRav0joUQGz1Wh1GU5Ms+fHD7ltVhFK35JaV4//ZNq8MoSpETWNa/NRaUlOHWz8NWh1G0Su9wWB1CSphEEBERERFRWphEEBERERFRWmyeQzCJICIiIiKyG7ZEEBERERFRemyeRMyzOgAiIiIiIiosbIkgIiIiIrIZdmciIiIiIqK02DyHYBJBRERERGQ3bIkgIiIiIqK0MIkgIiIiIqK0MIkgIiIiIqK02DyHYBJBRERERGQ3bIkgIiIiIqK0MIkgIiIiIqK0FHUSIYTo0f4MSCn35HJZuXC0pR2hYAh3K3fjj7b/YUrlf/3G3+DN4fMAgC9u+wKm/O/imZ37ca/nEwCAp5r3wOF05G8lCljrwWcRDAahKAq27vhySuXhUBgtzYdx7Hh71Liv9X4DY6NjM4ZTct2fO4hFZS5M/NSPvQOdM8r7njwOALh2M4gn/+SgPg0ABG6FTachcy3NrQgFgyivULB9x9aUys2G7dy2C+FwGOWKggOH9sE37sPO7bvh8SwFABxoOQAnj0NRDje3IhgMQakoxzaTujcrNxvmPTOI4XPDAICtO7bq9fxq72sYe2scHV08BiWTjW3hG/dh1/Y9WKrv8/u5z6dA/V4NQVHKE3zvRper37utOHa8TR/Pe2YQI0MjAICtX/kyz3syZPckYl6uZiyEWA6gTUq5BUCl9r5gvDl8HuVKOZ7rbIXL5cLl8beTlk/5pzB45iz2HtiNvQd26x+a32z4DTzX2YrnOlv5QUrR8NAIlAoFx463w+V2wTfuS1oeDoUx0P8tXPZdjho3HArDP+nPZ/hzxsOeOlz56SSaurfjvZsBPLC4Kqp806+vhXf8HJq6t+Pij8ew5sHVWPPgalz88ZieUMROQ+aGzw2jYrGCjq5jcJvt8yblZsNO9w3gQMsB9PR269MBQP3qenR0HUNH1zGeTMUYPjcCZbGCjq44xxuTcrNh/kk/zrz+BvYf2of9h/bp9RwKheG/ymNQKrK1LQCgfvXD6OhqR0dXO/f5FEx/r7Zp9Xg5abnZ965/0o83vuPFvoPPYN/BZ3jeMwtCZPbKl5wlEVLKS1LKCSGEG0CllPJSrpaVC28Once9VWrrQblyN84PX0ha/tdv/C2cLieOtrTjS5u36eO+638Xf3z86zMSEYpv5NwIllapV5DKlXKMDL2ZtNzhdGDdhsfhcEQfsE48/wLWbXg8P4HPMQ1Vq3Dpx+MAgImfTqK+qi6q/MKPxvDgR6vhWuDAgx+txqUfj2Plx+7DhR+NAQDO/+gtrPhYdd7jLkRDQyN6S4GilGNYu4qXqNxsWGPTGv2EqWKxgnA4DACYmvTjeOeJGSdlBAwPDcPj8QAAFEWZUfdm5WbDvG8MwuVy4nBzK7Zs+II+/YnOE3hi47r8rEyBy9a2AAD/5BS6uM+nbOb36kjScrPv3cE3zsLlcqL14LN4cuPMXhyUOiFERq98yVkSAeitEf0ABnO5nFy4Eb4Bp9MJAFjoXIhgMJi0PBgMolwpx94Du/FQ7Uq8OXweDudCAMDjGz6Hr3W9iCn/VH5XpECFw2G4XGr9OpwOBAPBtMojfOM+OBwOKBVKbgOeo9ylDly7qdZt4FYYd5W5o8q/d1VNMPq/0IXArTDe+f/8UYnDyo/dB3epM68xF6pwKAyX2wUAcDgcCMTu8yblyaYZPjeMmrpa/Qt+w6b16Oo8wZa5GGo9ascTh8nxxqTcbFgwEISyWMH+Q/tQu6oGw+dG1GOQk8egVGVrW0T2+fWb1qOr8wXu8ylQv1e140nc79345RHBQBBKhYJ9B59BbV3NjESQUlfUSYTWGtEAAEKI+thyIcRmIcQFIcSF117+Zi5DSdtCx0KEQiEAwI3QDf2Dk6jc5XJFtU68Pf4DvcuTw+nAQ7UrZ7RokDmHw4FgUK1f44lSquURJ57/Kvx+P1oPPovLvssY6PtWbgOfYwK3wlhUptatu9SB924Goso3/fpaXPzxGFY/vwHv3QxgzYOrMXDxDaxd8Wl0f+4gKj+sIHArZEHkhcf4pRwOh+GO3edNyhNN09Lciv0H1S41SoWid2OqXVXLL/UYaj1qx5OwyfHGpNxsmMvtiroi7vP5tBPYKRxubsVl32Wc7hvI45oVnmxtC3Wfb9f2eZ7IpkL9XtWOJ3G/d+OXR7jcrqgWi8sx3aIodUWbRMTcA/EegMrYcaSUJ6WUK6SUK9Zt/HyuQsnIQ3Ur8fb4DwAA42M+PTlIVH5v1Sf0YVP+d1Gu3B3VhWnKP4W7lbvztAaFrWZVjX7gGR8d1w9IqZZHdL/8dRw73o6tX/kylnqWYk3TZ3Ib+BzjHT+H5R9V72lY+bH79K5NEbGtDJGEY/XzG/R7IvrOfyf3gc4BdXU18Gn9ikdHx/VuSonK401zvPME1m9cp1/9NnbnmLzqh6KU53x9CkltXS18PrWOxkbH9EQgUbnZMI/How/z+/1QFAU9vS+io6sdW3dsxVLPUjQ2rcnjmhWebG0L4z7vv6puC0osW9+7S6uW6uNN+adQzuNNxor2nggAASFEjxCiH8BKKeXJHC4r6x6qXYlgMIijLe36ewB4ese+uOUP1a7ElH8KT+/Yhyn/FH5r9W/C4VyIp3fsw9GWdix0LNTnQ4nV1tUgGAii9eCz+nsA2LV9d8LyE50vwO+fwonOFxAOhS2IfG456xvCXWVu/WlLZ31DAKafyHTsO71oqFqFviePY+XHluGlv+uHa4ED3Z87iO7PHcSp868j+D63QypqV9UiEAiipblVfw+oT1qKV2427NXe1zAyNILDB1uxZcOT8J7xwuFwYOe2XWhpVltFI/MmVe0q9XhyWK9H9Xiyc9vuuOXxhvmv+rFz2274r/rR8MiMBnhKIlvbQt3nd+Owvs/XWLA2hSX+9+6ehOWx37u1dTXwT/qxa/se+Cf5OZgNu7dECCll3haWyPWf/dQegRShOz/0Ydz6+Q2rwyhapXcshNjoST4iZZ182YcPbt+yOoyiNb+kFO/fvml1GEVpQUkZALD+LbKgpAy3fs4LLFYpvcNh72enan7j1B9kdG78t4/9aV7Wjz82R0RERERkNzb/nQgmEURERERENmP3H5tjEkFEREREZDPzcpRDCCEaAQSg/o7bjHuWtYcjVQKAlPJ03PhyEx4REREREWUqFzdWawkEpJSD2nuzO9+f0pKHRUKIGU9XjWASQURERERkM/OEyOiVxEoAE9rfEwCMP8kAIcRmAOeFEJXaTzFMxM5Aj282K0dERERERNmXo0e8umPe3xXzfok27Jr2Uw2x4+uYRBARERER2cy8DF9CiM1CiAuG12bDbAMAFiVZ9BUpZQDARQCb443EG6uJiIiIiGwmha5JprSbpeP9yPN5TLdGVALwmpRHkgw31KTDPL6MoiMiIiIiopzJRXcm7YbpSu2GarfhBmuvodwdueHa7OlNEWyJICIiIiKymUxbIpKRUrZrfw4ahjUkKjfDJIKIiIiIyGb4Y3NERERERJQWu99zkDCJEEJcACCNg7T/pfa3lFKuzFFsRERERERFKVfdmbIlWUvEw3mJgoiIiIiIdHbvzpSwpURKGYy8oLY+tAHwau8XAdiUhxiJiIiIiIpKjn6xOmvSuSeiH8CTUBMJSCnfEUI0AOjIRiB3fujD2ZgNZaj0joVWh1DU5Ms+q0MoWvNLSq0OoagtKCmzOoSixvq3TukdDqtDIJuzdztEeknEXVriYHaPxKy9f/tmtmZFaVpQUoabPw9ZHUbRKrvDiVs/v2F1GEWp9I6FEI2VVodRtOTpCR77LRJJHj64fcviSIrT/JJS7vsWYvKcHekkEaeEEG8AuFMI8RmorRJ9uQmLiIiIiKh4FfqN1Top5TEhxCCAxwCsBLBHSvm9nEVGRERERFSk5kwSobkTwHtQb7KWScYlIiIiIqIM2P3pTCknEUKIFwFUQr3BWgB4WQjxhpTy6VwFR0RERERUjOZSS8RKKeUKw/uXtB+jYxJBRERERJRF9k4h0ksizqc4jIiIiIiIZqGgWyK0lgYJNRmqFEI0AbimFd8F4EpuwyMiIiIiKj4FnUQAeDgvURARERERka6gb6yWUgbzFQgREREREakKvSVCJ4RwAmiH+pjXiGtSyi9kPSoiIiIioiJm7xQCmJfGuGcBXIS6Tie194tyERQRERERUTGbJ0RGr3xJ5+lM16SUL2n9s96TUp4VQqzJUVxEREREREVrznRnwnSryiCAfiHEnwNYkv2QiIiIiIiKm91vrE65O5OUcrX2/zsANgMIAngwR3ERERERERWteRm+8iWdlgidlPISgEupjCuE2A31167XZrIsKxxubkUwGIJSUY5tO7amVB5vGAA4nA59mPfMIIbPDQMAtu7YCqfTkY9VKkitB59DKBhCuVKOrTu+lFL57u17EQ6HUa6UY99B9cfUB8+cxfDQCABg61e+BAfrPCWtB59FMBiEoijYuuPLKZWHQ2G0NB/GsePt+ni7tu9GOHwDilKOfQefyVv8c0X35lYscrgx8ZNJ7P3TNtNyAAjcCkeV7/q9zVj58U+iqfOP8hZrIcvWcR8AXu19DWNvjaOjqx2n+wbgPTMIAAgFg1j96Go8sWFdXtapkLQ0tyIUDKK8QsF2k/o3KzcbFgqF0XKgBR1dx6KmBdTvYrN5U27Pe9RtchgdXe0z5kuJFXRLhBDighDifJzXBSFEwl+sFkK4UWBdnobPjUBZrKCjqx0utwu+cV/ScrNh3jOD8FR7sP/QPgCAb9wH/6QfZ15/A/sP7cP+Q/uYQCQwMjQCpUJB+/GjWp1eTlo+0Pct7D/0DF58+Wv6OP5JP974jhf7Dj6NfQefZgKRomGtfo8dj/M5MCkPh8IY6P8WLvumt9VA37dw4NB+dL/8dX06St3D9/0qrvzkKpo6/wjvha/jgXuqo8rX1D6KixOjePKkepyJlLtKHVjykcV5j7dQZeu4D6gnTP6rfn3axqY16Ol9ET29L6KquooJhInhc8OoWKygo+sY3Kb1P7PcbFgoFMZA/wB8hmOQb9yH9RvX4cChfZia9M+YN+X2vCeyTYzfC5Q6u99YnazV42EA9XFekbJE2rRXwRgeGobH4wEAKIoy46THrNxs2NjoGKqXVQEAqpdVY2x0HN43BuFyOXG4uRVbNvDJuIkMn3sTnqp7AQCKUo6RoTeTlq9p+oyeJCgVCsKhGxh846/hdDnRevA5fGHjF/O7EgVs5NwIllYtBQCUm9S/WbnD6cC6DY/D4ZhO1GZuk3Ce1mBuaLj/U7g0MQYAmPjJJOo/+atR5Ss//klcuPIWAOD8D/8JK5bcBwBo+4M9aPuLnvwGW8CyddwHgBOdJ/DExnUzlhE5waKZhoZG4PGoxxNFKZ9R/2blZsOcTgee2LAu6hjkqfJAqVDU5M4/BU8Vt0GsXJ73TG+ThXlcI8qXhEmElDKY7BVvWiHEcgABKeVE1qPOoXAoDJfbCQBwOBwIBoJJy82GRT5AADA2OoZwOIxgIAhlsYL9h/ahdlUNhs/xqmw84XAYTpcLAOBwLpy5HZKUjwyNoKbuIbXOKxTsO/g0aupqMMIr4SkJh8NwubR92mnyOUhSbmZ4aAS1dTXZD3YOc5c6cO1GAAAQuBnCXY47o8qNicPKj38S7jInHrinGoFbYbzzb5P5DrdgZeu47xv3weF0QKlQZizjlZdeQWMTH2hoRq1L7XjucCBgWv/R5cmmMfKN+7Br+y7U8PhjKpfnPTQ7hd4SMRttACqFED0AlgshNseOIITYrHWLutD70is5DCV16glRCIB2oqQdpBKVmw1reKQe3jODONzcCv/kFBwOB1xuV1Tm7vOxWTUeh8OBUFA9kIVDN2ZuhwTlrQefwzPNT8Hh1Orc0GLhG387T2tQ2BwOB4JBbZ8OmXwOkpTHaj34LPY1sztZugK3wli00A0AcJc58V74elT5wPDrWFv32+je3IrKj1QgcDOEtj/Yg8qPVKB7cyuW31ONTfX/twWRF5ZsHfe7Ol+Af3IKh5tbcdl3Gaf7BgAA/kk/ypWZiQWpjBciwuEw3Kb1H12ebBojT5UHPb3dAKDfk0jTcnneQ7MjhMjolS85SyKklA3azdR7AFySUp40GeeklHKFlHLFhk3rcxVKWmrravWT+7HRMf2kP1F5vGl6el/U+wau/vRqeDwefTy/3w+FXypx1a56SD/hHxsd1xOBZOUnOr+KdRs+r18J9FTdq4/n909BUcrztQoFrWZVDS5r96GMj47rXZdSLTc60fkC1m143PTqLCXm/f53sbxSvc9h5cc/qXdtMlp9+PP6PRF9576N1Yc/j6bOP8KeP23DpXfG8NLgn+c15kKUreN+T++L6Ohqx9YdW7HUs1RvefCNX0b1fVV5XKPCUldXo9/HMDo6rndTSlSebJoIY/9+t9sFv38qF6tQ0HJ53kOzMw8io1f+4qMotatqEAwE9ScM1K5Smz93btsdt9xsWCgUxuHmVhxubsUjj66G0+lA7aoa+K/6sXPbbviv+tHwSLJbSopXTZ1ap60Hn9PfA+rTl+KVv9b7TYwMv4lnDx3BFzZ+EYNnzqKmrgb+ST92b98L/6Qf9Y88bM0KFZhavX6f1d8D6pOWEpWf6HwBfv8UTnS+gHAojNd6v4GR4RG0HnoOT278Q/0pNZSas2/9A+5y3Kk/gensW/8AAOjb8ccA1Buouze3ontzK079w7cRvMXuA5nI1nE/Hr/fz6uyCdSuqkUgENSfolS7qhYAsHPbrrjl8aY53nkCfr8fxztPIBQKw+FwoKW5FTu37cLYW2PsUmYil+c9ANDVeQJ+/xS6tG1CqbN7S4SQUqY2ohBOAO0AHpRSrhRC3ANgjZSyIxuBvH/7ZmqBUNYtKCnDzZ+HrA6jaJXd4cStn9+wOoyiVHrHQojGSqvDKFry9ATev33T6jCK0oKSMgDAB7dvWRxJcZpfUsp930ILSsrs/exUzVNDT2d0bnyk7rm8rF86LRH9UO9zeAfQf3SuIRdBEREREREVM5Hhv3xJ58fm7pJSviOEMGZFBZHJEREREREVErv/2Fw6ScQpIcQbAO4UQnwGwBYAfbkJi4iIiIioeOXzca2ZSDmJkFIeE0KcBdAE4CEAe6WU38tZZERERERERUrY/PlH6bREQEp5CcClHMVCRERERESYQy0RQogfAjDeD1EJ9fcfVmY9KiIiIiKiIjZn7omQUn7c+F4I4QawN9sBEREREREVu3w+aSkTaXVnMpJSBoQQD2QzGCIiIiIimlvdmS5gujuTAOAC748gIiIiIsq6OdOdCcDDsQOklMEsxkJERERERADmzaGnM52UUj6Ws0iIiIiIiAiA/Vsi0klxLgghduQsEiIiIiIiKgjptEQ0AKgXQjwN4BqAIADJR7wSEREREWWX3VsiEiYRQoiNUsqXtbdr8xAPEREREVHRm1fgj3hdC+BlgDdRExERERHlS0G3ROTTgpIyq0MoamV3OK0OoaiV3rHQ6hCKljw9YXUIRY3HfmvNLym1OoSixX2fkin034lYIoTYmWgEKWVHNgIJ/zsbOqzi+AUXPrh9y+owitb8klLWv0VY99aaX1IK8XsfszqMoiT/8kcAgPdv37Q2kCK1oKSMxx4LFUryPBd+sfqXch4FERERERHp5onC/p2ICSnl3rxEQkREREREAOx/T4S9UxwiIiIioiIkMvyXdL5CNAoh6oUQm5OM15aoPFkS4U0aCRERERERZdU8ITJ6JSKEaAQAKeWg9r4+znj1ACoTxpeoUEp5LGEkRERERESUdTlqiVgJIPJYwgkAy2csV4hKwzhxsTsTEREREZHNZNoSIYTYLIS4YHgZuy25YxZzl8miK6WUSZMI2/xOBBERERERqUSGT2eSUp4EcDJOcQDAovjLFPWRrk7JMIkgIiIiIrKZHP1OxHlMt0ZUYub9z9e0+yHcACqFEMullJfMZsTuTERERERENpOLG6ullKehJgf1ANyGG6y9WvklbdgizOz6FIUtEURERERENpOr34mQUrZrfw4ahjXEjJOoSxQAJhFERERERLYzLzfdmbKGSQQRERERkc3wF6uJiIiIiGhOYUsEEREREZHNZPqI13zJWRIhhFgOoB9A5LFQm6SUgVwtLxeePXQEoWAQ5Uo5vvyVL6VU/uyhIwAAh2OhPmzPV/YiHLqB8opyPNP8VP5WoMC1NLeq9VuhYPuOrSmVmw3buW0XwuEwyhUFBw7ty+s6FJpc1jm3Q3LZqv9QKIyWAy3o6DqmT2s2jBLr/sJzWORwYeJfJ7H3m0dNywEgcDOEvd88igcqq9G/++u4NDEKANj0tacQvBnKa8yF6nBzK4LBEJSKcmwz2ffNys2G7dy2GwDgcjmxXzvOHG5uBQA4nA7TeRerbB1v4n0HAIDT5dKP9y2G7WC2PJrJ7vdE5DrFOS2lXKu9AjleVlaNDL0JpaIcbc8fhcvtwuXxy0nLL49fxuPrP49nmp/ClH8Kl8cv41v9/w37Dj6Dr7/8x/p0lNzwuWFULFbQ0XUMbrcLvnFf0nKzYaf7BnCg5QB6erv16chcLuuc2yG5bNV/KBTGQP8AfL7pY5bZMErs4ft/FVf+9cdoav8i3gtfxwOV1VHla1Y9iotX3sKTLz4NAHr56XOvo6n9i2hq/yITiBQNnxuBslhBR1c7XKb7/sxys2Gn+wZQu6oGHV3t8FR74D0zCO+ZQXiqPXpCETvvYpWt402874DaVbXo6DqGqmoPvGe88I37sH7jOhw4tA9Tk35uhxQJ9den037lS66TiEohRJvWKlFQRoZGsNSzFABQrpRjZPjNpOVLq5ZCqShHOBTGlP9dLK1ais+s/c9wOB0AoJaFw/ldkQI1NDQCj1a/ilKO4aGRpOVmwxqb1sCp1X/FYoX1n0Au65zbIbls1b/T6cATG9bB4XDo05oNo8Qa7v81vUVh4idXUX//p6LKV378flz44T8BAM7/8PtY8fFPAgAqf6UCRz+/d0bSQfENDw3D4/EAABRFmbHvm5WbDateVgXfmJpI+8Z88FQtxdjoGKqXVQEAqpdVY2x0PI9rZl/ZOt6YDateVoVxbTuMj/ngqfLAU+WBUqEgFArD75+Cp8qTx7UtXCLDf/mSyyQioP1/BECbEKIyh8vKunD4BlwuJwDA4XAgGAimVH55/DL2fOUpPFS7csY8R869iZrah3Ic+dwQDoXhcrsAqPUbiK1/k/Jk0wyfG0ZNXW0eoi9M+apzbgdzuah/ypy7zIlr4QAAtbvSXY47o8qNicPKj98Pd5kTAa3l4cjA19H2+F7c85GKvMZcqNT9OMH3rUm52bDIiemu7bvhcDqgVChRicPY6BgvYGiydbwxGza9HXbp2wFQW4F2bd+FmrqavKzjXFC0LRFSyglDNyYvgPrYcYQQm4UQF4QQF159+bVchZIRh2MhgkH1CyEcnv6QJCtfWrXUtOvSs4eO4Onmp/RWCUrM4Zz+IgmHw3DH1r9JeaJpWppbsf/gPv1qOM2Ujzrndogv2/VPsxO4GcIihxuAmlC8F74eVT5w7nWs/dXfRvcXnkPlr1QgcDOEd34yqXdj8v7jd2e0XpA5dT9O8H1rUm427HTfADzVHvT0vgiX2wXvmUE0PFIP75lBHG5uhX9yiq1xmmwdb8yGne4bQFW1Bz293XC7XfCe8QIAPFUedmlN0zyIjF75iy9HYrowLQEwETuOlPKklHKFlHLFExvX5SqUjNTU1eCy1n94fHRc77qUqNx434TL7cKUfwoA8MLzX8Xj6z8PpaI8T9EXvrq6Gr3/9ujouN5cmqg83jTHO09g/cZ1+tUQMpfrOud2SCyb9U+z5/3+32N55TIAaktDpGuT0ermz+n3RPT9w7ejujAt+ZWPYuInV/MTbIGrrauFz6f2kR8bHdO7KSUqNxsW28oQDKontz29L+r3RKz+9OqcrkuhyNbxxmzYzO0QiroHwu12wa+dH1FiQszL6JW3+KSUuZmx2n2pDcA1AAEp5Z5E44f/PZibQGbhhee/inD4xownLbU9f9S03D85hW+88k2EtANX2/NH8Y1XvgnvdwbhdKrNrp9p+s+oX/2wNSsUh+MXXPjg9i2rw5jheOcJhEPhqCc57Ny2S3+6jFl57LBXe1/Dmdff0LueNTatQcMjDeYLtMj8klLb1H+u6tzvn7LldrBT3QPZqf/IsIH+AaxZuwYbNq2H0+kwHWa1+SWlEL/3MavDiOvo5/fq3ZQiT2fq2/01NLV/Ea4yJ9o+vxeAmnAMnHsd93ykAm2P78W1cDBqGjuSf/kjAMD7t29aG4imy7AfG5+01NHVHrc83jQRHV3tCIXCONF5AgBQu6oWDY/M6BRhiQUlZZYfe7J5vDGbT0RH1zH4J/145eXX9PMjq58SN7+k1N6PPdL0XfnTjM6Nm5b8QV7WL2dJRLrsmEQUC7smEcXCbieyxYR1by27JxFzmd2SiGJjhySimBVKEtE/8f9kdG68tvKzeVk//tgcEREREZHN5PNJS5lgEkFEREREZDP5fNJSJphEEBERERHZjN1/sZpJBBERERGRzbAlgoiIiIiI0iJy+pvQs8ckgoiIiIjIZtgSQUREREREabH705ns3U5CRERERES2w5YIIiIiIiKbmcfuTERERERElA67d2diEkFEREREZDO8sZqIiIiIiNLCR7wSEREREVFa2BJBRERERERpmcd7IoiIiIiIKB1siUiR4xdcVodQ1OaXlFodQlFj/VuHdW8t+Zc/sjqEoragpMzqEIoWjz2UDJ/OlKLQv1+3OoSi5fyFO/HB7VtWh1G05peUsv4twrq31vySUrx/+6bVYRSlSPIgHlticSTFSZ66gls/D1sdRtEqvcNhdQgpYUsEERERERGlhU9nIiIiIiKitPAXq4mIiIiIKC28J4KIiIiIiNLCeyKIiIiIiCgtbIkgIiIiIqK0sCWCiIiIiIjSMs/mT2eyd3RERERERGQ7bIkgIiIiIrIZdmciIiIiIqK08MZqIiIiIiJKC1siiIiIiIgoLWyJICIiIiKitDCJICIiIiKi9BRzdyYhRCOABu3tHillIJfLy7bnDh1FKBhCuXI3vvSVP0q5/Juv/AnGR304+vxzCYfRTC3NrQgFgyivULB9x9aUys2G7dy2C+FwGOWKggOH9sUdRtGyVf+hUBgtB1rQ0XVMn5b1by5bdZ7K5+B03wC8Z7wAgGAwhEceXY0nNqzLz4ra3OHmVgSDISgV5dhmsh3MymOH+cZ92LV9D5Z6lgIADrTsh9Pp0D4Ph9HR1Z7XdSp03ZtasWihCxM/mcTeP5tZd92bWgEAgZshvdxsGKWu9eCz6j6tlGPrji+nVB4OhdHS3Ipjx9sAAAN934L3jUEAQCgYwupPN2DdhsfztxJziN1bInL2OxFCiEoAj0kpt2ivQK6WlQsjQ29CqSjH0eefg8vtwuXxyymVh0Nh+CenosY1G0YzDZ8bRsViBR1dx+B2u+Ab9yUtNxt2um8AB1oOoKe3W5/ObBhFy1b9h0JhDPQPwOeb/syw/s1lq85T/Rw0Nq1BT283enq7UVVdxQRCM3xuBMpiBR1d7XCZboeZ5fGmqV/9MDq62tHR1a4nEAP9A7jsu2y2aIrj4ftW4cpPrqLp+Jfw3o0AHrinOqp8Tc2ncXFiFE++pF6QeOCeajxwTzXa/rIHT760D5UfqZgxDSU2PDQCpULBseNt2j59OWl5OBTGQP+3ovbvNU2fQffLX0f3y19HVbWHCcQsCCEyeuVLLn9srhHANSFEjxDCm8Pl5MSbQ2/iXs+9AIBypRxvDp9Pqfyrx7+Gz6//XNS4ZsNopqGhEXi0K3iKUo7hoZGk5WbDGpvWwOl0AAAqFisIh8Omwyhaturf6XTgiQ3r4HA49GlZ/+ayVeepfg4ivGe8qKr25Hz9CsXw0DA8HrU+FEWZsR3MyuNN45+cQlfnCT2pmP48LMzX6swJDfd9CpfeGQUATPzkKurv+9Wo8pVLPokLE28BAM5f+SesqLwP33tnDO/82yRcpQ5U/vJifO+dsbzHXchGzo1gaZV6HClXyjES8zkwK3c4HVi34fGo432E98wgPFU8zsyGyPBfvuQyibgLwBUp5RYAXiFEfQ6XlXXh8A24XE4AwELHQgQDwaTll8cvw+FYCKWiXB/PbBiZC4fCcLldAACHw4FAbJ2blCebZvjcMGrqapMOo9zUvxnW/7Rs1Xm6n4Pel15FY9OanK1XoVHrTz2eOxyOmcd7k3KzYZETqfWb1qOr8wX4J/15XIu5xV3mxLUb6nYI3ArhroXuqPJI4gCoCYW7TN0WD9xTjf7tf4zBt/4hr/HOBeFwGC6XdhxxmnwOkpTHevXl17Cm6TO5CbZIFHMS8R6AS9rfEwCWx44ghNgshLgghLjw6suv5TCU9DkcCxEMhgAAN8I39C/oROVfff5rmPK/i+cOHcXl8bfxrf6/MB1G5owHpXA4DHdsnZuUJ5qmpbkV+w/u06/GxhtGqmzXvxnWf7Rs1Xk6nwP/pB+KwosaRmr9qcfzcDg883hvUm42TKlQ9G5MtatqZrRoUOoCN0NYtFDdDu5SJ967EYgqHxj5DtbWPoruTa2o/EgFAjfVbfG9d8aw+lm1+8zD963Ka8yFzuFwIBjUjiMhk89BknIj/6Qf5TzOzFquujMJIRqFEPVCiM0mZW4hxHJtnLZE88llEnEJ04lDJdREIoqU8qSUcoWUcsUTG9flMJT0PVT3EN72vQ0AGB8d17suJSr/2stfxdHnn8OXtn8RS6vuxWfW/r7pMDJXV1ej96MfHR3Xu2ckKo83zfHOE1i/cR2UCkWf3mwYTctm/Zth/c+UrTpP53PgG/eh+j72FTeqrauFz6d2PxobHdO7KSUqNxtmvJfCf9UPReG+ninvW9/F8nuWAVBbGiJdm4xWP/u4fk9E39C3o+6BeO9GAJW/vDg/wc4RNatq9Ps7x0fH9a5LqZYb+cYvo3pZVe6CLRK5aInQHnoEKeWg9j62p1ATgBVSytNa+YxEIyJnSYQW3BIhRD+AJZFgCkVN3UMIBoJ47tBR/T0A7P3K0wnLKXO1q2oRCATR0tyqvwfUJ8zEKzcb9mrvaxgZGsHhg63YsuFJeM94TYdRtGzVP6CevPr9fhzvPIFQKMz6jyNbdZ7q5wAA/P4p0/7Lxax2VQ2CgSAO6/VXAwDYuW133HKzYQ6HAzu37cbh5lY4tNYIAOjqPAG/X71XIhTi/UCpOPvWOdy10K0/bensW+cAAH3bvwoAcJU60L2pFd2bWnHq3LcRvBVG4GYI3Zta0bf9q1i55D68dPbPLYu/ENXWqft068Fn9fcAsGv7noTlJzpfgN8/hROdLyCs7d9TPM5kRY66M63E9IX9GT2FtAv8J7W3po0AenxSygxXLbtC/37dHoEUIecv3IkPbt+yOoyiNb+klPVvEda9teaXlOL92zetDqMoLSgpAwCIx5ZYHElxkqeu4NbPmVBapfQOh72fnaoZvX4po3Pj+xY9uAWAsQVBTwyEED0AeqSUl7RWiAYp5Z7YeWhPWd2j3dtsij82R0RERERkM5neJK0lDCfjFAcALEphNo2JEgggt/dEEBERERFRBnJ0Y/V5AG7t70oAM/oXCyEapZTt2t8zHowUwSSCiIiIiKgIaPcoV2pdmdyGG6y92v/1ANqEEBeFEBeRoNWC3ZmIiIiIiGwmV7/5EGllADBoGNag/T8IIKWbpZhEEBERERHZTD5/OC4TTCKIiIiIiGwmlR+OsxKTCCIiIiIim2FLBBERERERpYVJBBERERERpYXdmYiIiIiIKE1MIoiIiIiIKA1siSAiIiIiorTwnggiIiIiIkoLkwgiIiIiIkoLuzMREREREVFa7N4SIaSUVscQYZtAiIiIiGjOsvfZuebqjSsZnRsvXrgkL+tnm5aI92/ftDqEorWgpAwf3L5ldRhFa35JKevfIqx7a7H+rTO/pBQAv3utsqCkDKJBsTqMoiW9fqtDSAm7MxERERERUVrs3p2JSQQRERERkc3YvSVintUBEBERERFRYWFLBBERERGRzbA7ExERERERpYlJBBERERERpcHeKQSTCCIiIiIi27H7jdVMIoiIiIiIbIdJBBERERERpcHeKQSTCCIiIiIiG7J3GsEkgoiIiIjIZux+TwR/bI6IiIiIiNLClggiIiIiIpsp2h+bE0JsBrBWe7sIwCkpZXuulpdNh5tbEQyGoFSUY9uOrSmVxw473TcA75lBAEAoGMTqR1fjiQ3rsHPbboTDYShKOfYf2pfX9SoELc2tCAWDKK9QsN2k7s3KzYaFQmG0HGhBR9cxfVqzYRSN9Z9/2arzeMMAwOF0JBxW7HZu2wUAcLpcOBBzXFaP5V4AQDAYwiPasTyR2G3hG/dh5/bd8HiWAgAOtByA0+nI/ooUsGx870a82vsaxt4aR0dXe9zvYkque+tRLHK6MfEvV7H35edMywEgcDOkl5sNo8zZPYnIWXcmKeVJKWWDlLIBwIVCSSCGz41AWaygo6sdLrcLvnFf0nKzYY1Na9DT+yJ6el9EVXUVntiwDqf7BnCgZT96el/U50XThs8No2Kxgo6uY3Cb1v3McrNhoVAYA/0D8Pku69OaDaNorP/8y1admw3zjfuwfuM6HDi0D1OT/rjDit3pvgHUrqpFR9cxVFV79IQhQj2Wd6Ont1s/licSb5vWr65HR9cxdHQdYwIRI1vfu4B6rPFf9evTmn0XU3IPP/ApXPmXH6Pp8JN4L3QdD3x8WVT5ml/7HVz857fw5Im9AIAHPr7MdBjNbTm/J0II0QjgYq6Xky3DQ8PweDwAAEVRMDw0krQ80TTeM4PwVKtljU1r9C8PZbGCcDic8/UpJENDI/qVOkUpn1H3ZuVmw5xOB57YsA4Ox/QXtdkwisb6z79s1bnZME+VB0qFop5U+afgqfKYDit21cuqMD6mJr/jY764deI940VV9XTZzm27sHPbLpzuG4gaL942nZr043jnCSZuJrL5vXui8wSe2LhuxjKM38WUXMODv45L//wWAGDiX66ifvmvRZWvvPd+XPjB9wEA59/+PlZ84n7TYTQ7QoiMXvmSjxurn5JSnszDcrIiHArD5XYCABwOB4KBYNLyRNO88tIraGxaM2M5w+dGUFNXk6vVKEhqPboAqPUYMK376PJk01DqWP/5l606jzcf37gPu7bvijrWmA0rZpGkYdf2XXA4HVAqFNPxel96VT+Wv9r7Gh55dDU6uo5hfCw6KTDbFpHkecOm9ejqPAH/pB80LVvfu75xX9xtGO+7mMy5y5y4Fg4AAAI3grjLeWdUuTFJWHnv/XAvdJoOo7ktp0mEEKISwEQul5FtDqcDwUAIABAOT38ZJCqPN41/0o9yZebB7HBzK/YffIZN2jHUelS/PMLhMNymdR9dnmwaSh3rP/+yVefx5uOp8qCntxuA2s0m3rBidrpvAFXVap243a4Z3ZkA9ViuKOX6+8mrfgydG8HxzhNwOB043nkCxztP4HTfgOm2UCoUvRtT7araGVfai122vne7Ol+Af3IKh5tbcdl3WW8livddTPEFboawyOEGALgXuvBe6HpU+cDffxtrf/130L31KCr/w2IEboRMh9HsiAz/5UuuWyKWAzgfr1AIsVkIcUEIcaH3pVdyHEpqautq4fOpV5bGRsf05tJE5fGm8Y1fRvV9VVHTd2lNrfGudhWzuroavc/86Oi43iUgUXmyaSh1rP/8y1admw0zdptxu13w+6dMhxW72G6lweDMEx/fuA/V91Xr76uqPahbVYPtO7ZGvRqb1iTdFpNXoxMSyt73bk/vi+joasfWHVux1LNUb3kw+y6mxLwX/w7L/4/7AKitCpGuTUar9/4X/f6Hvv/5V3GH0WyIDF/5keskohJAIF6hdvP1Cinlig2b1uc4lNTUrqpBMBDEYe0JJrWr1Cb/ndt2xy2PN43f74/qA/5q72sYGRrB4YPPYsuGL+hPjCBV7apaBAJB/ekxtatqAUw/OcWsPN40xztPwO9X+yCHQuG4w2ga6z//slXnZsMcDgdamluxc9sujL01hsamNabDit0TG9Zh+NzwjDqJbAMA8Punoo7ljU1rMHRuRK9Lo3jbYue2XWhpboVDa42gadn83jUT+11MyZ393ndxl/NO/WlLZ7/3XQBA3361FdNV5kT31qPo3noUp/72rxC8GTIdRrNj7xQCEFLKPC4uvvdv37RHIEVoQUkZPrh9y+owitb8klLWv0VY99Zi/VtnfkkpAOD92zctjqQ4LSgpg2hgjwSrSK/f3s9O1Vz/2U8zOje+80Mfzsv68cfmiIiIiIhsJp9PWsoEkwgiIiIiItthEkFERERERGmwdwrBJIKIiIiIyIbsnUYwiSAiIiIishm73xORj1+sJiIiIiKiOYQtEURERERENpPPX5/OBJMIIiIiIiLbYRJBRERERERpsHcKwSSCiIiIiMh27H5jNZMIIiIiIiLbYRJBRERERERpsHcKwSSCiIiIiMiG7J1GMIkgIiIiIrIZu98TwR+bIyIiIiKitAgppdUxRNgmECIiIiKas+x9iV/zwe1bGZ0bzy8pzcv62SmJKGhCiM1SypNWx1GMWPfWYv1bi/VvHda9tVj/1mL9E7szZc9mqwMoYqx7a7H+rcX6tw7r3lqsf2ux/osckwgiIiIiIkoLkwgiIiIiIkoLk4jsYb9A67DurcX6txbr3zqse2ux/q3F+i9yvLGaiIiIiIjSwpYIIiIiIiJKC5MIIiIiIkqZEOJDVsdA1mMSMQtCiCVWx1CshGqFEGKnEMJhdTzFRqv/aiHEk0IIp9XxFBshxD1CiCNWx1GMtH3/QSHEcSHEAqvjKUZCiOVCiGYhRKnVsRQTbd9/SAjx3wD8X1bHQ9ZjEpEhIYQA0CeE+I+G95Q/XwXwuwB+CcCXhBCLLI6n2LQCeBhAJYAtFsdSjJoA/Ebk+EP5oV2wOAfgPwFokVK+b3FIRUcI8VsAVkopDwG4bXU8xUIIUQlgAEA9gKsAfmxtRGQHTCIydz+AcQBfBgDJO9TzRrvyfVZKeQBAJ4AyAAutjap4CCHmAfgfUsoXAPwJgF+zOKSiodU9ALwNoBvA5ywMpxj9LwCnAPwLgDohxCohRInFMRWbKgCLhRDNAL4phPiw1QEViQCAjVLK5wB8G8Cd1oZDdsAkInPfl1J+DsBdQojlAFsj8igM4LtCCCGl/CmAxVATCcoDKeX/BvA9IcR8AJ8FEBBC/JoQ4hctDm3Ok1L+by2R+E9Sym8AeF8IsUMIwSQ6D6SUPwPwD1C7cvwigCcA7BBCuCwNrLjcBuAA8LcA+gF8TgjxS5ZGVASklNeklNe0t58H8G8Az3uKHZOIzEU+OKehtUZQfkjVT6WUUgihAJiUUvqsjquYSClvQT1+/HcAJ6B273jU0qCKhJbETQkh9gH4KNR6v2VtVEXlHwH8HYB/AtAGwAlgqZUBFZlLAEqgtsa9ATWpWGFpREXC0BJ6EUAjwF4YxY5JRIa0L3JIKf8YwBIhxO8CKLc2qqL0ywD+qxCiTgixnlfD80dKeUtK+V0p5UWo3TyCVsdUDIQQdwD4jwDeBbAJwE8ALLM0qCIipfx3KeUxKeWElPKHAD4M4N+tjquIXAQwCfW+iBCARQCuJZ6EsiFy3gO1Ne5nfKgG3WF1AIVMy8rnA3BBvRL7P62NqCg9AeBBAOcBtEkp/5fF8RQNIUQ5gCUAHgNwA8D3rI2oOEgpfy6EeFhK+XMAEEK8BMBvcVhFRdv37wPw+1DvjXvb0oCKiLb//1cATUKIzwL4ZwBsic6vDwFwA/gPAELWhkJW4i9Wz5L2iL9flFLyKmyeaVdkmwB8x9BXk/JICFEBICylDFgdS7ERQswzXBmkPBNCfAxAiMcea2itzvO11gjKI+1+uBIp5U2rYyFrMYkgIiIiIqK08J4IIiIiIiJKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIiIiIiJKC5MIIqIUCSGuCyG8Qoh+7f+2WczrivZ/ZeTvbEk0PyGEWwhxMcn0ScfJZFwiIpo7mEQQEaVnrZRyrZSyAQCEELtnMzPtR8uWJBpHCFE5m2UQERFlG5MIIqLMnQLQkIfl9OdhGURERCljEkFElLnHAHi1Lj1eIcRuIUQ/oLZQCCEuasPd2rAebVgPgEXasKjuQEKIzdo4F7V5tAFYHpl/OvNOhTYPrzat21BUGTPPyPgzlk1ERMXnDqsDICIqMP1CiMjfl6SU7drJdD0Ar5RyrRBiOYAGKeWDWlekNi25qJRSPggAQoim2Blr020xjOOWUgaEEPWG7lMZzTsew3w3A2gDsMVQtkUr8woh6gFci102gD2pLouIiOYOJhFEROlZK6UMmAwPSCnbtb8fA7Ao0ioBoBJAANHdkq6ZzOMxqF2kAABxlpPpvE1pyUA9ZnbLMs6jH8ByAHeZLJuIiIoQkwgiouyIPXE/IqU8HXkzmyc5mcjKvLVWjZcAbAIwAUMrRBrLdmeybCIiKmy8J4KIKDsChr9PwXBCrl3t9wJYq713w/y+hVNQWxqM08EwzWzmbWYFgFNSykuY2apgnMcWAINxlk1EREWILRFERFkmpbwU+S0JbZBXu3dirXYT9QWoV/7Npusx3Gg9CPWeg0EAZ4UQF6SUWzKZdxx9AC4KIRqgJkHG1pQJ7YbqFQAGtUQDscsGcDLFZRER0RwipJRWx0BERERERAWE3ZmIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgtTCKIiIiIiCgt/z+Td/NezMfiNAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Build the plot\n",
    "plt.figure(figsize=(15,5))\n",
    "sns.heatmap(matrix, annot=True, annot_kws={'size':10},\n",
    "            cmap=plt.cm.Greens, linewidths=0.2)\n",
    "\n",
    "# Add labels to the plot\n",
    "class_names = ['1', '2', '3', '4', '5', '6', '7']\n",
    "tick_marks = np.arange(len(class_names))\n",
    "tick_marks2 = tick_marks + 0.5\n",
    "plt.xticks(tick_marks+ 0.5, class_names, rotation=25, fontsize=10)\n",
    "plt.yticks(tick_marks2, class_names, rotation=0, fontsize=10)\n",
    "plt.xlabel('Predicted label', fontsize=12)\n",
    "plt.ylabel('True label', fontsize=12)\n",
    "plt.title('Confusion Matrix for Random Forest Model - Kano River\\nTraining data - Kano\\&Yoshii River, 75\\% training and 25\\% testing, Accuracy score is 96.09\\%')\n",
    "plt.savefig('data/kano_river/out_img/predicted/confusion-matrix-kano&yoshii-data.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5d4c612a-7f5e-486a-898f-220a567e7803",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the trained Classifier to the prediction dataframe\n",
    "preds_yoshii = rfc_model.predict(Yoshii_classified[features])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ae636ed0-ef8b-4d9d-bbf5-9e15164ff7e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our classification accuracy is: 97.14533764790879%\n"
     ]
    }
   ],
   "source": [
    "# View accuracy classification (cross-validation) score\n",
    "print('Our classification accuracy is: {cv}%'.format(cv=accuracy_score(Yoshii_classified['label'], preds_yoshii)* 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2089e347-2cb5-4295-acf5-e15f52ac4117",
   "metadata": {},
   "source": [
    "*Visualizing confusion matrix*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3d9cd274-00ae-421f-a33c-3c924dbc97b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get and reshape confusion matrix data\n",
    "Matrix = confusion_matrix(Yoshii_classified['label'], preds_yoshii)\n",
    "matrix = Matrix.astype('float') / Matrix.sum(axis=1)[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "681bdafe-9c86-4961-b00d-08c0e562f909",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxEAAAFZCAYAAADn86JOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABTs0lEQVR4nO3de3gU150n/O9BfrIgqO42THZnoxaeyLnQAidjGaxLkt2ZtYQnkzf7zA4XPzu5mIsBJ54EMAhsBwskZIO4GMm5GEGI7Zl937wglJn3zSQeLLGTnSRIMoJMAlLjxCgxaiV5kzX0BbAnG/a8f1RVq7pVfammq6ta/f3w9IP6nO6qU6dOVdevzqkqIaUEERERERFRtmY4XQAiIiIiIiouDCKIiIiIiMgSBhFERERERGQJgwgiIiIiIrKEQQQREREREVnCIIKIiIiIiCxhEEHkQkKIGiFEnxDinPb/ttucTp+F71QJIS7nMr8U07sshOgxSe8QQtzWPaZzKWsudWIyjWvaNHq0/7tznZbJtPNW94ZpXjO0Jf3ly/d8sihH3tpCunoSQviEEOdyLGOPEKIjRfryLL5vOm9jW83UbtO1r3xvn0REuRJ8TgSRuwghagD0AGiSUo5paVX63xan1QdgQy7fzRftgCcspbwvKf0cgCop5Z0Zvp/TsqeZ3m3XiRDiGoD3SinD2vv1AO6TUm7IQ/kuSynvvt3pJE0zobx2yGY93W5bSJ5WqnrSAqTTyfPJcro+AD+Huj717a8RwHYpZVOW389p3oZp2Na+iIjyhT0RRO5zFEkHubdxwJvXA/Db0G88i6sFSsNZfnfKmevbZEed9ANYnOdpFpts19PttAXbaQfuewAYe5e6ATh5AM/2RUSuwyCCyH2qpJT9ZhlCiPWG4SiNWppPG6rSrf3foaVvA1ClfbYmeZhF0tCKc4bhLsmfy3qeaRwH8JDh/UNIOug0zOOcPtRGm64+/GibNl/97x5jWbWyrNf+rkkxpCShTsyWL3keGZZLP/PcAcNBZ4plSVlnhrRuAHMN6enqvsPw/zZ9fpnKm1T2VNNPWH7tb2P7SG4zCespw2yzaQtTypWhnhLKZ6UOzEgp9wGYK4RYri3PSSnlWIr6SqgLbRI+k+3R2FazHm6V3L6SptOjt2Ptvb5NJ68vS22aiCgrUkq++OLLJS8AVQAup8irAdBjeH8OgE97SQA+Lf2a4TOXDX/7AJxLzoN6gLLc8Jn453KZp0m5L5uU5XKq7wFYD6DbOM+kZZAAtqVYpnNamc/pZUtVnjTL90fGeaSYxjUAfdpLAqhJ8bn4sqSqMwCNAPqM07ZY99Kw/vrMypJU3j4A27KY/jZDOfoM7bM7uc0kr6fbaQtpypWunpLL58umPBnKWqPV2+UM5TLbfszWc7xMmcqXrn0lTWc5gA5j+dLUR9o2zRdffPFl9XUHiMg1pHq2c26K7A1Qz+Lq+qEeWPUDGJOT492vWpxtH4BuIcQSqAccxu/nc54nhTqMZUz7fgIhRJU27UzjzsNSPVNsZgXUg7t1Mrvx/2bL9ycZ5hGfl5QyrJ3ZrQJwXs9IsyxmddaExDPxenq2dW+szzEYztCblddQxu400zcu/0NQz8rrZawCsB2JbSYMa9K1hVTLvQTm9WRWvtsmpTwvhBg2zDNVucy2n9vZHnUp25ehjCe1no7tUOtB7+Uxq49s2jQRUdY4nInIfcaMQzgy0A8Yw7nOTKpDp/SD3WyGWOQ6T31cudnwFf0s6jASx6KbSXdQ5tPyb+dAUp9GtrZDPRsNIOOyhG+jXECKus8yYLIy/eTl3yOlXKG97suhzSRL2RYylCuVhPKl+6A2zEgfkpSpnYSRvi3MTVEX4QzTtSKhfZkY05ZjuZTypJZmVh+5BjNERKYYRBC5zzqoZzZ9eoLhwNR4cad+1tgKnzY9H7QDbaHdVUdKuR3qga/P8Pl8zBNA/OLwuQAa5dRrPhYDOC6lPA+TACBpnHs4zWyOArgPwBLjWPE0zJbvexnmkUBbrn6hXY+BDMtiog9qD4q+nPoBc97qPoV00w8b0o8bP6ddA5DQZozXfWQz4wxtIVW5UtXTlPJlmPcGKWWT9rJygb1puTJsP7fNpH0l64YaaOg9FanqI5zPchERcTgTkctowyhWAOgxDG06LqXcJxIvGN6jDX/yZTndsBDivHbx53lMHnQ0avMD1OElYcN3+m9nnib08dnJTgA4J4RowtSzv/0ATmtDS7anmrA2fGO7tpzrtO88kO4svdnyAfhF9osT1wH1IPdIhmVJVYYVWhmGoa4DO+rebL4Zp6+1xx7DRcN9AMLGNqPVeXw9yexuRWraFlKVC+oZd7N6MivfkawqwYI09bU+1faTR8b2lawfaoCzQitnQeqDiIjPiSAiIiIiIks4nImIiIiIiCxhEEFERERERJYwiKCciMmHW/UJIS5rL/19qgsAk6dRpT8c6XY+kws7ppnlfK8lXzCt1Z0v9bdyntd67YFTUx7+pY2ZnnLHFy19eXK6yedMH5ZlXF+Z1p1WF336+G3tlqNTpuMWYvLhYTVJ6Vm1dwvzcWy7sCqfZShU/SZNuyD1qG3jmR7IWPLc0q7NONE+idyOQQTlREq5T7/DCYCTUB+opd/xJKuL+LQ7mtx9u5+xSxa3f8zH9HsANOXxFp1GTdrtHR8yyVsHYL1xGYV6W1mf4TaRlhnXV5brTr8NZRPUi5G7LXw3Z4YAS39d04MnIYQ0BMT6E5troN4NaAWAJw3T2Qb1QupU87Hchty+XWRLTD5R+5xIfOq0o/Vr/E4h6lGb3xjUB8NRGlbXh0h8grjxCfSm23bSd33aZ81OsviS9gk5tU+i6Y5BBFFq2dy/Pidaz0Mf1INoK7eZtDJ9X6oeDi1o2YPE5xjo9+53Sj/U26PaTkp5RHvewX0AHgAwbAiezhsCYv2uO1VQ70I0hslb4/qgPsArnGZWtrUhN9OChsta/a5AYj04Xb+FXicboG1nyWexKXdacKafgFgB9fbOmbZtow7te2aehHoHMF2u7ZNoWmMQQXmnncXp04Y86Wca9bNF54ThnvL6kBj9b0OXcUe2n9HyOgzz6DbrYjZ8rxuGh1elKFsHgBp9OVJ97jachtp7k/AU2nT1lGK5jWfi4md7tR+2Kqg/rNsNn/cZPrMP6pNtl2vLeNJwy8qEaWrDDPQze/qtI33ZrK9sKkMrVwcmD7aM0+kxHnyJyeFS2/TyaJ+f0u6ylFBHKZyH2h6Mz1JYn67XLbkNuXi7sNrmTLejJGN63Rhuz5qu1+C261dLS24TCe02xTrJSz2moT8E7iSSegVF4lnzbWZpyduRof1n1Z5STLNHJPYOTRlCZLbNpylv8v7CrGwJ68asopLWh9k+J6FeoZ6I0fd3V03aWMptW7sN8ZT9hDaNGiQ+Idxy+yQqCVJKvvi6rRfUg79thvc+ANKYZshbD/XgWf/cuaTv+LT31yx8pgpAj6EsjSbzbQTQZ3h/LV3ZtPfnUixvwucs1tU1qD98Pammn6KezJa7Rl9uvbz6ZwzzMq6XmuRya2nXoJ41TjlNrV6XG8qT7fpKt4x6XfRp06pJakP6dJYD6DCWT/u/z7D+9WcOmLa7NGWITyepXD3asjca0pdr6T7te3p91KSZ/rmkZXLVdpFDm8u4HaWY9uUC1K9Zm0hotynWSV7rMal8VYYy1RjrS3ufXBazNF9S2uWk8mZqT2bTXG5YpoRt3vC55G0+VXnN9hcJZTNbNynqy7g+pqw7k2XsSJq3sT1N2bZTTGNbUpq+f4nPP9f2yRdf0/3Fh82RXcJSPdMNIH52pxFAU5rvjMnJruFUD+hK9Rmf9v9bhr+NmpB41in+3SzLlvXnsjAmpdygnQ1MOJuVZh5my70B6tNpdf3ad09qZ+6OQD3zqa+HJ6EOYYqT6oOphjFZN6mm2Qf1KdpLoB6cXU1RJqtWSPVBZT1QDy7OJ39ASnlSOyu8XVue49r/c8Vkj0P8qbzGdpeF+FATgyNSffowhBDXANyplwPq2WQIIR6SUm7X5n9VCBHWv5OB27YLq20u5XaUilZHxmEjdtWvWZvYjsR2G07xXV3O9ZjCBkyeLT+vnWmvkWoPpN6WoeWHhRBmaenml017mjJNqPuIDrN8g+Rt3qxsqfYX/UllS7W9ppM8/3BSfj+A7UKIPVCHQSYPFTPbttMSWo+ntq4Seo3ytP0TTSsczkR2MR6k62ePh5F+px7OYrpTPiPV4RJztQPnJdLChcHZli2bz2nDHPRu/XQ/kvqPzToAHWJyCEm6eYQzLoxqrjbU4rz2ozasDSNYDuCqTBo+ZZh2uoPBuVLKfkwemOhDK7ItUza2Qz3zl4o+HGa5Yf3ukepF2SukOv4ZsB7MrMTk8AQAQNLBwFUx9W4sjQCOa+utSqrDIhqRHVdtF3lqcykJdcjTHmO7s7l+E9pEinabTjg54Xb2L1DPXi/RThj0QL3A2uxGB7nKpT3p+rW6bjRbphzqTqcPcUveFs2215QyzV9bL9uhDlnST0AYr2OYsm1n4UkAVdq6Xg7gSZF0UfZttk+iaYVBBNklbPh7MYDj2oFE3u94pB1cHpeJF2om64N2NlTb+es/dGnLZjgLmHEZpJQb5OQFoxkvlpaTFzcfzXYeSXqQeCG0cbyuTy+T9pkO7e+cpimEqJLqnVO2Q/2h9mUxraxp9dUvUo8174Z6wKAfjB43ltMQtIVzmHf8O0KIRn1aejsxCbyatLRU1wMkSDqTbCyfG7YLq2VItR2Zzb8bSdf92Fy/U9pEcrs1BOy+5OmkmX429Zjqe2HDgfMKqHWnt3G9Ry1hPiZpgLa96QeuhtmEDX+nWpepptmtlcd0X2Wyzb9iMp10+yBj2VJtrymlWndGUsqTWr0OI7EXSc8PJ38nHUPw2aQtxx6TAMtS+ySazjiciQrhBNTbdzYh81nvXFyFekZf/5E6mdy1LKXsF0Ks0C7aG8bkD2e6svUDOK0N99luxzJIKfcJ9R7yNRnKYvbdfqE+Z0I/S7dHOxg/op351NNPQr0IWn+f8o5QqaYp1Isn9QOoMeS3F0LXAfUg1exixX6oBywrtHKe15ZRv+DS9HtCvWD0PrODCTF5602jYQBHDQfIDyR9Zz20s7xavfi06aRaV8ltyMjx7cJqGdJsRwm0eloJYLEQQk9eBxvr1zBE0NgmwsZ2qw3BSbdOzKSsx3TtCybDabR2e1Uf0qT1XurbWr82RMYs7by2XOdhMuRPY7ouzeYDYLuW3pimDhqT6u57Kcpmtr/wmSx38rrJNOwwef7h5A8YpjdmDPBSbNsJ60sLcvULwZdkEyDm0D6JpjUhpXS6DES3Rahjey9LKY9oP15HoZ6Ry/l5BzQ9CCG6s+yByXZ6jdowC/19DbRrTVIMFXNMMW4XbqzfdPWY7/ZVaEKIPu2sO2XBje2TyEnsiaDpwIfJs25hIcQY2NVMqrw+Jdh4AKG9P4/U95p3mg9Ftl24tH59SF2PRfsUai0gCjtcjKLi0vZJ5Bj2RFDRM5wd9GlJY8V8dpAoH7hd5Md0rEehPuOhCdqwJqfLQ0TFiUEEERERERFZwrszERERERGRJQwiiIiIiIjIEjddWM1xVURERERkN5H5I84TTf6cjo1lX6ggy+eaIEJ88i6ni1Cy5LffxDu3bjpdjJI1s6wcb9+64XQxStKsstls+w5i23fOrLLZAMD6d8isstmsewfp7Z9uj2uCCCIiIiIi0gh3d5gwiCAiIiIichuXX7nMIIKIiIiIyG3YE0FERERERJa4O4ZgEEFERERE5DrsiSAiIiIiIkt4TQQREREREVnCnggiIiIiIrLE3TEEgwgiIiIiIteZ4e4ogkEEEREREZHbuDuGYBBBREREROQ6vCaCiIiIiIgscXcM4fabRxERERERkdvY2hMhhPABOCqlXGHnfArh8GPPYu4cH8Z+fQVPvLzXNB8AwtejeOLlvVj34H/Fio98AgAwV/Hh+Pe/jf3f6i5omYtR2852RCMRVFT6sXnLxqzyU6UBgOJR4mnRaAxtLW040Lm/QEtTfHbvbEckEoW/sgKbTOrfLD857eSJXvSd6gcARCMRLP34UqxeuwpbN20DAHi9HjzduqNwC1WE8rUdsM1nLx9tX/fisZcwcmEUBzr3ITgaRPPm7VgQWAAAaGl7Gh6PUpiFKiJ27nv6TvVj8MwgAGDjlo2sfxN21j+QuE2QBS6/sNq2nggtgFgPoMaueRTKAx/+CC7/6k2s7Pg83opdw713L0rIX9bw5zj3xgU8+tWnAAD33r0IR099E0tbPo2lLZ/G8Bs/YQCRhcEzg6ic78eBzv3w+bwIjgYz5pulBUeDWPPIKrS07sDEeAjB0SCi0Rh6e3oRDF5yaOncb/DMEPzz/TjQuQ9e0/qfmm+WtnzlMnQfewHdx15A9cJqrF67CidP9KKuoRYHOvchsDAQ/6GhqfK1HbDNZy9fbR9QA7fQlVDC9xuXPoADnftwoHMfD2BN2LnvCY2HcOqVV/F06w483bqD9W/CzvoHzLcJypLI8VUgtgURUsqwlHIfgLBd8yiUpj/+GM5fvggAGPv1OBo//JGE/CXv/xCG37gAADj7s59g8fvuiefpAQZlNjAwhIB2ts7vr8DgwFDGfLO0QHUA/kq/uuMKTSBQHYDHo2D12lVQFP6ApDI4MIhAIAAA8Pv9U+rfLD/dd/pO9SOwUM1buKgawRH1wDY4EkSgekEhFqko5Ws7YJvPXj7bftfBLqx+ZFXC90PjE+g82DXl4IxUdu57+l7th9frwe6d7diw9nOFWJyiY2f9A+bbBGVJiNxeBcJrIrLgm+PB1esRAED4RgTzlDsT8o2Bw5L3fwi+2Z543pMrPo+jp75ZuMIWsVg0Bq/PCwBQFAXhcCRjfqrvqEMImlFbX1vAJShual2qbVdRFERM6z8xP913vnH0G1i+chkAIFCt/qA0b94GxaPAX+m3fXmKVT63A8pOvtp+cDQ4pX3rQdyadWvQefB5hMZ5RjaZnfueSDgC/3w/nm7dgbqGWgyeSTxAJnvr32ybIAtmiNxehSpeweZUxMLXo5g7R/2B9s324q3YtYT83jPfxYqPfAKHH3sWVX84H+EbUQDAe/9dJcZ+faXg5S1WimdyRxSLxeDTDorS5af6TqA6gO5jhwEgPhaW0lPrUm27sdjkQWm6/FTfCY2HUOGf/NE4eaIXgYUBdB97AV6fl8OZ0sjndkDZyVfbV4OECeze2Y5LwUs4eaIX/kp/fBhTXUPtlLO8ZO++x+vzJpwxDwbZG5TMzvo32ybIglIdzpQNIcR6IcSwEGIYb153sihp9f3L91GjXQex5P0fig9tMlra8un4NREnfvAdAEDN3ffg7M9+XLiCFrn6+tr4+O2LF0fjwzPS5ZulGYcM+HxehEITBVqC4lZXXxf/gR25OBL/4U2Xn+o7wdFLWHhPdfy7sVgsYVqRCM+Up5Kv7YCyl6+2333sBRzo3IeNWzZiQWABlq9clrA/Cl0Jwe/nGdlkdu57AoFA/HOhEOvfjJ31b7ZNkAWlPJxJCNEBoEoI0aFdaJ1ASnlESrlYSrkYd82xsyi35fSPf4h5yp3xOzCd/vEPAQAntn8NAOCd7cHhx57F4ceexfHv/wMiWk9E1R9WxnslKLO6hjqEw5H4nZXqGuoAAFs3NafMN0tTFAVtO9uxdVMzRi6MxHdahw52IRQK4dDBLkSjseTZl7y6hlpEwhHsjtelOhRMv6uSWX6q74RCoYSx+KvXrsLgmSFs3bQNIxdG+UOSRr62A4BtPlv5bPvJFEXB1k3bsHtnOxStN4IS2bnvqWuoRehKCFs3bUPoSghNDzYWbLmKhZ31T7fJ5T0RQkpZuLmlIT55lzsKUoLkt9/EO7duOl2MkjWzrBxv37rhdDFK0qyy2Wz7DmLbd86sstkAwPp3yKyy2ax7B80qm+3ue6dqxMMfzOnYWL78ekGWj0+sJiIiIiJyG5eHOgwiiIiIiIjcpoDXN+SCQQQRERERkdu4/B6qDCKIiIiIiNyGPRFERERERGSJu2MIBhFERERERK7DnggiIiIiIrKE10QQEREREZElLu+JcHmMQ0REREREbsOeCCIiIiIit3F3RwSDCCIiIiIi15nh7iiCQQQRERERkdu4/JoIBhFERERERG7j7hiCQQQRERERkdsI9kQQEREREZEVDCKyJL/9ptNFKGkzy8qdLkJJm1U22+kilCy2fWex7TuL9e8c1j1l4vIYwj1BBBERERERqWa4PIpwTRAR+d1bThehZHnfNQ/iMx9wuhglS/7tT3Hj9zGni1GSZt+h4J1bN50uRsmaWVbO+neI3gPH+ncG276ziqUHmsOZiIiIiIjIEgYRRERERERkCYMIIiIiIiKyxOUxBIMIIiIiIiK3YU8EERERERFZwiCCiIiIiIgsEWAQQUREREREFri9J2KG0wUgIiIiIqLiwp4IIiIiIiKXcXlHBIMIIiIiIiK3mWFTFCGEWA4gDKBKSnnEan68fLaUjoiIiIiIciaEyOmVYZrLAUBK2a+9b0zKbwQwpuWPCSFqUk2LQQQRERERkcvYEUQAWAJgTPt7DEBykDAMoEcLHqqklOdTTYhBBBERERGRywiR60usF0IMG17rDZP1Jc1mnvGNlDIMoBtAD4D70pXP1msihBA9UAs7JqXcYOe87LCntQORSBQV/vfgC48/ljH/0ujreHLLl/DBwAcBAF/a9QQUj4JYNIZndu3F3ueeKfQiTAuHV7dh7hwvxn4zjieOHzDNB4DwzWg8/8QXuuAr92DsN+N49MWWgpZ3umjf9QyikSgq/BXYuOWLWeXHojHs3tmOfYc6Ej4HAIqimE6HJrXtbEc0EkFFpR+bt2zMKt8sbeumZsRiMVT4/Whp3YHgaBBbN29DILAAANDS1gKPRyncghWBfNV92852AIDiUeJp0WgMbS1tONC5v0BLU3zsrH+zNEpkZ/0n748oe7ne4lW7jiHVtQxhAHPTzLMRQL+Ucp8QokMIsVxKedLss7b1RGhRzzopZZOhUEXjtYGzqKiswN7nnoHX58Wl0dezyv9PTX+Kvc89g73PPRMPIP7u5P+D14Ovm82GMnhgYQMu/+YKVn55I966Hsa9d1Un5C9b8iDO/WIkHijce1c11v3pQ1h3bAeWdqyOT4OsGRoYgr/Sj32HOuD1eREcvZQxPxaNobfnW7gUnPxs/6l+BKoD2LHrSwAwZTo0afDMICrn+3Ggcz98Pi+Co8GM+WZpJ0/0oqWtBd3HDse/BwCNSxtxoHM/DnTuZwCRJF91HxwNYs0jq9DSugMT4yEER4OIRmPo7elFMMi2n4qd9W+WRonsrP9U+yPKjk3Dmc5isjeiCkBfUn6NYQjTHqQJOGwLIqSUR7QuEQC4jKndJ642NPAaFmg9ChX+9+C1wbNZ5U+EfokvP/fVeFCheBR8ds2noSj80c5F06IGnP/FKABg7DfjaFyUGBAsqboHw2MXAABnxy5gcdU9OPpPxxG5GQMAXP7NFfjKWfdWDZ4ZQqBaPWvt91dgaGAoY77iUbBq7cOYY2jrIxdHUb1IDfwWLqrG6MXRAi1B8RkYGIr3FPj9FRhMqnOzfLO05SuXxYOEyvl+xGLqtjAxHsKhg108iDKRr7oPVAfgr/QjGo0hFJpAoDoAj0fB6rWr+BuQhp31b5ZGieys/1T7I8qOHUGE1qtQpZ3c9xkusNaDiSPacKhGACvdcHemJgD9BZpXXlyPXYfH6wGgDsOIhCMZ8xVlDgBg1SOfxVcOfQ0T4xOFLfQ05Cv34Or1MAAgfDOGeXN8Cfl64ACoAUVywNC06CPoHzlTiKJOK7FYDB6vF4AaCCe3/0z5OmPgMHJxlD8gacSiMXh9Wp0qCsLJdW6Sn+k7g2cGUVtfFz+AXbtuDToPdiE0HrJ7cYpKPus+OBpE8+Zm1NbXFnAJipvd9c91kl6h2r++P6Ls2dQTASnlPillv5RynyGtSfs/rHUE9KcLIIACBBFCiG4AGwy9Esa8+IUfL339ZbuLYskcZQ6ikSgA9YBJ31jS5evDmxSPgvvrl0zpvSDrwjejmKsFDr5yBW9pAYWu9+wprLj/z3B4dRuq/m0lwjcnD1IPr27Dhm+0xHslKHuKoiAaUX8UjD8W2ebrGh9sRP+r/Wjf9QwmQhM8G5uGMRiLxWLwJde5SX6677TtbMfTu3bA41Hgr/THhzHVNdRNOdNY6vJZ94HqAIduWGR3/XOdpFeI9m/cH1H2cr2wulBsDSKEEB0AOqSUY2b5WqSzWEq5eNUjD9tZFMtq6+/HJe06htGLwfjQpXT5xusmJsYnUOGvKFyBp6m+i2dQ80fqcJglVffEhzYZLe1YHb8m4sTQdwEAex/aio5/OIqf/3a8cIWdRuoaauPXL4xcHI0PXco23+iFr38tfk3E0j9rsqnExa++vjY+bv7ixdH4UIF0+am+c+hgF9Y8sgr+Sj8AJAxhGr8Sgp/7pgT5qntjPft8XoRC7I3Ohp31z3WSmd3tP3l/RNmzqyciX+y8sHobgEYA3UKIPv3hFsXi/voliIQj2NPaEX8PAE88/qWU+YoyB088/iXsae3AHGVO/Dtffu6r8WslYlGeFbfi9MgZzJvji9+B6bQ2NOnEF7oAAN5yBYdXt+Hw6jYcH/wuIjdjaP7EI2hc2IDuNW14dfuLWLbkQcfKX6xq62sRCUfid1bSu6a3bd6eNr/r4POYCE2g6+DziEVjiEVjaN/1DNp3PYOlf9YEhWehUqprqEM4HInf3aSuQe3237qpOWW+WdqLx17C0MAQdu9qx4a1j6LvVB8URcHWTc1o29kOReuNoEn5qntFUdC2sx1bNzVj5MIIlq9cBkA9iAqF1GtSovwNmMLO+k+1TmiSnfVvtj+i7Lk9iBBSyoLNLJ3I795yR0FKkPdd8yA+8wGni1Gy5N/+FDd+zwMLJ8y+Q8E7t246XYySNbOsnPXvkJll5QDA+ncI276zZpaVF3DQT+4qn/3TnI6Nx5/6p4Isn63PiSAiIiIiIusKeX1DLhhEEBERERG5TCGHJuWCQQQRERERkcsIuDuIKNRzIoiIiIiIaJpgTwQRERERkctwOBMREREREVnCIIKIiIiIiCxxeQzBIIKIiIiIyG3YE0FERERERJYwiCAiIiIiIksYRBARERERkSUujyEYRBARERERuQ17IoiIiIiIyBIGEUREREREZAmDCCIiIiIissTlMYR7ggjvu+Y5XYSSJv/2p04XoaTNvkNxuggla2ZZudNFKGmsf2ex/p3DuqdM2BORpcjvrjpdhJLlfddcvH3rhtPFKFmzymZDPPxBp4tRkuTLr7PtO2hW2Wy8c+um08UoSfoBLOvfGTPLynHz9zGni1GyyovkxB2DCCIiIiIisoRBBBERERERWeLyGIJBBBERERGR27AngoiIiIiIrHF5EDHD6QIQEREREVFxYU8EEREREZHLcDgTERERERFZ4vIYgkEEEREREZHbsCeCiIiIiIgsYRBBRERERESWMIggIiIiIiJLXB5DMIggIiIiInIb9kQQEREREZElDCKIiIiIiMiSkg4ihBDd2p9hKeV2O+dlhz2texGJRFHhr8AXHn8s6/y/+cZ/w+jFUex97lkAwBOPPwUA8Ho9eHLnE4UpfJHavbMdkUgU/soKbNqyMat8s7S+U/0YPDMIANi4ZSM8HgXRaAxtLbtxoHNf4RZoGji8qhVzZ/sw9ttxPHHigGk+AIRvxuL5Jx7rgq9cwdhvx/HoSzsLWt5ilo/2HxwNonnzdiwILAAAtLQ9jYlQaEqax6MUbsGKQNvOdkQjEVRU+rHZpO7N8lOlAYDiUeJpWzc1IxaLocLvR0vrjgItUXGxs/7N0shc+65n1P2JvwIbt3wxq/xYNIa2ne3Yf6gj/rnmzeohn9frwY5dXypM4achtwcRM+yasBCiBkCHlHIDgCrtfdF4beA1VFT6sfe5Z+H1eXFp9PWs8mPRGCbGQ/HP/V3P36O2/n7sfe5ZLKhegNOv/veCLkcxGTwzBP98Pw507oPX50VwNJgx3ywtNB7CqVdexdOtO/B06454ANHb04tLwUsOLV1xeqC6Hpd/M46VX92It66Hce9d1Qn5yxY/iHO/GIkHCvfeVY11f7IS617cgaX718SnQZnlq/0DQOPSB3Cgcx8OdO6LBwtmaaQaPDOIyvl+HOjcD59p3U/NN0sLjgax5pFVaGndgYnxEIKjQZw80YuWthZ0HzscnxYlsrP+zdLI3ODAEPyVfuw/1KHtTy5lzI9FY+jt+VbCb2vviW+hrr4W+w91IFAdQN+p/kIvyrQhRG6vQrEtiJBSnpdSjgkhfACqpJTn7ZqXHYYGzmJB4AMAgAr/e/Da4Nms8r9y6Kv47JrPxD8XWBjAJW1DuzR6CQsCHyzQEhSfwYFBBAIBAIDf78fgwFDGfLO0vlf74fV6sHtnOzas/RwAwONRsHrtKijKnAIuUfFrWvQRnP/FCABg7DfjaFzYkJC/pOoeDP/8IgDg7NgFLH7vIhz93glEbsYAAJd/Mw5fuaewhS5S+Wr/ABAan0Dnwa6EAyazNFINDAwhoPXS+P0VU+reLN8sLVAdgL/Sj2g0hlBoAoHqAJavXBYP2irn+xGLxQq4ZMXBzvo3SyNzQ2eGsKBardMKfwWGktaDWb7iUbBq7cNQlMkTE9WLqhEcDSIWjSE4GkRA+w5ZJ4TI6VUotgURQLw3ogdA0YWh12MxeLxeAICizEEkHMmYf2n0dcxRFFRUVsQ/t6BaDRqe3PKlKXmUKBaNwetTDzgVRZlS52b5ZmmRcAT++X483boDdQ21GDyTuCOk7PnKFVy9oa6H8M0o5s3xJeTrgQOgBhTJAUPTwgb0j54pSFmLXb7av/5jvmbdGnQefB6h8ZBpGk1S61HfnysIm9Z9Yn6q76jDyZpRW187ZT6DZwZRW19n56IUJbvrP906oUmxWAxe/bjGY7IPypCv04OG5se3Q1EU+Cv9NpZ6eivpIELrjWgCACFEY3K+EGK9EGJYCDH80tdftrMols1RFEQj6gYSi12P76zS5X/l0FcxEZrAnta9eD34U/xdz9/j73r+HguqF+ArR5+H1+flcKY01J1SFIC2s0qqc7N8szSvz5twdjYY5JnXXIVvxjB3troefOUevHU9nJDfO3wKK5b8GQ6vakXVuysRvhmN5x1e1YoNL7XEeyUovXy1f3+lPz5kqa6hNj4EITmNJhkPiGKxGHymdZ+Yn+o7geqA6dCltp3teHrXDg4lM2F3/adaJ5RIURRE9OOaqMk+KEO+rvfEtxCoDuDw178Gr8/L4Uy3oWSDiKRrIN4CUJX8GSnlESnlYinl4lWPPGxXUXJSW78El4I/BQCMXgzGhy6ly//K0eex97ln8debH8MHAx/Af1nxF4jFrid8LxqJgszV1dfFD/hHLo7EA4F0+WZpgUAgnhYKheD38yxIrvou/hA1f7QQgNrToA9tMlq6f038mogTr70CANi7cis6vnMUP/8tz3hnK1/tP2EI0xW1/Zul0aT6+loEtTHdFy+OxofJpMs3SzPWs8/nRSg0AQA4dLALax5ZxTOyKdhZ/6nWCU1V21CLS9p1EKMXR+NDl7LN1yUP2eNxT+5K9poIAGEhRLcQogfAEinlERvnlXf319+PSDiCPa174++ByTstpcpP9tk1n8bQwGt44vGnMHpxFP9lxV/YX/giVddQi0g4gt3anTTqGtSu562btqXMT5UWuhLC1k3bELoSQtODaidY58EuhELquPBolGfHs3F6dADz5vjid2A6PToAQL37EgB4yxUcXtWKw6tacXzou4jcjKH5zx9B48J6dK9qw6vN38CyxQ86Vv5ikq/2rygKtm7aht0726FoPQ9maTSprqEO4XAkfhefugZ1yNHWTc0p883SFEVB2852bN3UjJELI1i+chlePPYShgaGsHtXOzasfRR9p/ocWEJ3s7P+zdLIXF29uj9p3/VM/D0weaelVPldB59HKDSBroPPIxaNYdXahzE4MITmzdsxcnEUy1b+pQNLMz24vSdCSCkLNrN0Ir+76o6ClCDvu+bi7Vs3nC5GyZpVNhviYV5w7wT58uts+w6aVTYb79y66XQxStLMsnIAYP07ZGZZOW7+nieznFJ+h+Lue6dq/uT4p3M6Nv7eQ/+tIMvHh80REREREbmNy58TwSCCiIiIiMhl3P6wOQYRREREREQuM8PdMQSDCCIiIiIit2FPBBERERERWTLDpiBCCLEcQBhAldndU7XHNFQBgJTyZMry2VI6IiIiIiLKmR23eNUCCEgp+7X3Ux4GDeBJLXiYK4SY8pw3HYMIIiIiIiKXmZHjK4MlAMa0v8cAGB8ODSHEegBnhRBV2kOhx5InYCwfERERERG5yAwhcnoJIdYLIYYNr/WGyfqSZjMv6f3dWtpV7aHRyZ+P4zURREREREQuk+uF1dp1DlOuddCEAczNMInLUsqwEOIcgPUA9pl9iEEEEREREZHL2HRh9VlM9kZUAegzydeDDB/UoMMUhzMREREREbmMHRdWaxdMV2kXVPsMF1j3GfJ9+gXXZndv0rEngoiIiIjIZew60y+l1Icn9RvSmtLlm0kbRAghhgFIY5I+fe1vKaVckmWZiYiIiIgoC3Y9JyJfMvVEPFCQUhARERERUZzbn1idtqdEShnRX1B7HzoA9Gnv5wJYV4AyEhERERGVlFxv8VooVq6J6AHwKNRAAlLKnwshmgAcyEdBvO/KdLcpstOsstlOF6GkyZdfd7oIJYtt31kzy8qdLkJJY/07p/wOxekikMu5ux/CWhAxTwsczK6RuG1v37qRr0mRRbPKZuOdWzedLkbJmllWzvbvkFllsyE+eZfTxShZ8ttvct/jED14YP07g/t9Z/HkUX5YCSKOCyFeBXCnEOIvofZKnLCnWEREREREpavYL6yOk1LuF0L0A3gIwBIA26WUP7KtZEREREREJWraBBGaOwG8BfUia5nhs0RERERElAO3350p6yBCCPEC1Mdj90C9FuLrQohXpZRP2VU4IiIiIqJSNJ16IpZIKRcb3h/VHkbHIIKIiIiIKI/cHUJYCyLOZplGRERERES3oah7IrSeBgk1GKoSQqwEcFXLngfgsr3FIyIiIiIqPUUdRAB4oCClICIiIiKiuKK+sFpKGSlUQYiIiIiISFXsPRFxQggPgH1Qb/Oquyql/FzeS0VEREREVMLcHUIAMyx89jSAc1CX6Yj2fq4dhSIiIiIiKmUzhMjpVShW7s50VUp5VBuf9ZaU8rQQYplN5SIiIiIiKlnTZjgTJntV+gH0CCH+bwB3579IRERERESlze0XVmc9nElKuVT7/+cA1gOIALjPpnIREREREZWsGTm+CsVKT0SclPI8gPPZfFYIsQ3q065X5DIvJ+3e2Y5IJAp/ZQU2bdmYVX40GkNby24c6NwX/5xZGplr29mOaCSCiko/NpvUuVl+qjQAUDxKPE1dD2040Lm/QEtT3HJp/2ZpWzdtAwB4vR483bqjcAswTRx+7FnMnePD2K+v4ImX95rmA0D4ehRPvLwX6x78r1jxkU8AAOYqPhz//rex/1vdBS1zsbNzP0Tm8lXnZvt57vszy9f+vu9UPwbPDAIANm7ZCI9HMU2j7BR1T4QQYlgIcTbFa1gIkfaJ1UIIH4p0yNPgmSH45/txoHMfvD4vgqPBjPnRaAy9Pb24FLwU/5xZGpkbPDOIyvl+HOjcD59pnU/NN0sLjgax5pFVaGndgYnxUMK6CXI9ZCWX9m+WdvJEL+oaanGgcx8CCwPoO9Xv0BIVpwc+/BFc/tWbWNnxebwVu4Z7716UkL+s4c9x7o0LePSrTwEA7r17EY6e+iaWtnwaS1s+jeE3fsIAwiI790NkLl91braf574/s3zt70PjIZx65VU83boDT7fugMejmKZR9tx+YXWmXo8HADSmeOl56XRor6IzODCIQCAAAPD7/RgcGMqY7/EoWL12FRRlTvxzZmlkbmBgCIHAAgCA318xpc7N8s3SAtUB+Cv9iEZjCIUmEKgOGNYDd2DZyKX9m6UtXFSN4Ij64x4cCSJQvaCwC1Lkmv74Yzh/+SIAYOzX42j88EcS8pe8/0MYfuMCAODsz36Cxe+7J56nBxhkjZ37ITKXrzo3289z359Zvvb3fa/2w+v1YPfOdmxYq9793yyNpo+0QYSUMpLpleq7QogaAGEp5VjeS10AsWgMXp8HAKAoCiLhiKV8sk6tUy8AtU7DpnWemJ/qO8HRIJo3N6O2vraASzB95NL+zdL0A6fmzdugeBT4K/0FXIri55vjwdXrat2Hb0QwT7kzId8YOCx5/4fgm+2J5z254vM4euqbhSvsNMH9UOHls87Junzt7yPhCPzz/Xi6dQfqGmoxeGbINI2yV+w9EbejA0CVEKIbQI0QYn3yB4QQ67VhUcPHjn7DxqJYp3gURMJRAEAsNrmzyjafrFPrVN15xWIx+EzrPDE/1XcC1QF0HzsMAPGxmJS9XNq/WdrJE70ILAyg+9gL8Pq8HM5kUfh6FHPnqHXvm+3FW7FrCfm9Z76LFR/5BA4/9iyq/nA+wjfU+n/vv6vE2K+vFLy80wH3Q4WXzzon6/K1v/f6vAm9E8Fg0DSNsieEyOlVKLYFEVLKJu1i6u0Azkspj5h85oiUcrGUcvHadWvsKkpO6urr4o195OJIfCPINp+sq6+vjY9bvXhxNN5VnS7fLM04ntPn8yIUmijQEkwfubR/s7RYLJbwvUiEZwut6PuX76NGuw5iyfs/FB/aZLS05dPxayJO/OA7AICau+/B2Z/9uHAFnUa4Hyq8fNU55SZf+/tAIBBPC4VC8Pv9pmmUvRkQOb0KVz4yVddQi0g4gt3a3TXqGtTuaP1OM6nyOw92IRSaQOfBLkSjsZRpNFVdQx3C4Uj8jiZ1DXUAgK2bmlPmm6UpioK2ne3YuqkZIxdGsHyl+kzEQwe7EAqFcIjrIaNc2r9Z2uq1qzB4ZghbN23DyIXR+Lqg7Jz+8Q8xT7kzfgem0z/+IQDgxPavAQC8sz04/NizOPzYszj+/X9AROuJqPrDynivBFlj936IpspXnQPm+3nu+9PL1/6+rqEWoSshbN20DaErITQ92GiaRtlze0+EkFJm90EhPAD2AbhPSrlECPFeAMuklAfyUZC3b93IriCUd7PKZuOdWzedLkbJmllWjrdv3XC6GCVpVtlsiE/e5XQxSpb89pvc9zhkZlk5ALD+HcL9vrNmlc12971TNU8OPJXTsfGe+mcLsnxWeiJ6oF7n8HMg/tC5JjsKRURERERUykSO/wrFysPm5kkpfy6EMEZFRRHJEREREREVE7c/bM5KEHFcCPEqgDuFEH8JYAOAE/YUi4iIiIiodBXydq25yDqIkFLuF0KcBrASwP0AnpBS/si2khERERERlSjh8vsfWemJgJTyPIDzNpWFiIiIiIgwjXoihBBvADBeD1EF9fkPS/JeKiIiIiKiEjZtromQUr7P+F4I4QPwRL4LRERERERU6gp5p6VcWBrOZCSlDAsh7s1nYYiIiIiIaHoNZxrG5HAmAcALXh9BRERERJR302Y4E4AHkhOklJE8loWIiIiIiADMmEZ3ZzoipXzItpIQEREREREA9/dEWAlxhoUQW2wrCRERERERFQUrPRFNABqFEE8BuAogAkDyFq9ERERERPnl9p6ItEGEEOIRKeXXtbcrClAeIiIiIqKSN6PIb/G6AsDXAV5ETURERERUKEXdE1FIs8pmO12EkjazrNzpIpQ0tn/nyG+/6XQRShr3Pc5i/TuH+33KpNifE3G3EGJrug9IKQ/koyA3f389H5OhHJTfMQfv3LrpdDFK1syycta/Q1j3zppZVg7x8flOF6MkyVeuAADbv0O473FWsQTP0+GJ1X9geymIiIiIiChuhiju50SMSSmfKEhJiIiIiIgIgPuviXB3iENEREREVIJEjv8yTleI5UKIRiHE+gyf60iXnymI6MtYEiIiIiIiyqsZQuT0SkcIsRwApJT92vvGFJ9rBFCVtnzpMqWU+9OWhIiIiIiI8s6mnoglAMa0v8cA1EyZrxBVhs+k5JpbvBIRERERkcqmW7z6kt7PM/lMlZSyP9M1GQwiiIiIiIhcRuR4dybtWgfj9Q5HpJRHtL/DAOam+W6jPtQpEwYRREREREQuk+tzIrSA4UiK7LOY7I2owtTrn69q10P4AFQJIWqklOfNJsS7MxERERERuYwdF1ZLKU9CDQ4aAfgMF1j3afnntbS5mDr0KQF7IoiIiIiIXMau50RIKfdpf/Yb0pqSPpOuNwMAgwgiIiIiIteZkeNwpkJhEEFERERE5DJ8YjUREREREU0r7IkgIiIiInKZXG/xWii2BRFCiBoAPQD020Ktk1KG7ZqfHdp3PYNIJAK/34+NW76YVX4sGkPbzt3Yf2hfwmdfOvYyRi6OTEknVdvOdkQjEVRU+rF5y8as8rNJC44GsXXzNgQCCwAALW0tePUfX0XfKfWOZpFIFA9+fClWr11VmAV1uXyth2g0hraWNhzo3J/wXQBQPIrptEuVXW0fALZuagYAeLxetLTuAGC+bsjc4S/swVzFh7FfX8ET39hjmg8A4RvReP6Jp14AAFyNhfHol58sXGGLXL62g62bmhGLxVDh98fbvNl2QIX93Z0IhaakeTxKgZa0eLn9mgi7Q5yTUsoV2its87zyanBgCP5KP/Yf2gevz4vgaDBjfiwaQ2/Pt3ApeCnhs7FoDKHxUCGLX1QGzwyicr4fBzr3w2dW1yb52aYBQOPSRhzo3I8Dnfvh8ShYvnIZuo8dRvexw6heWM0AQpOv9RCNxtDb04ugYTsIjgax5pFVaGndgYnx0JRplyo72/7JE72oa6jDgc79qF4YQN+pPtN1Q+Ye+OOP4vKv3sTKZz+Ht6LXcO/dixLyl330z3HujQvxQOHeuxdh3cf/Cn0/+j5WPvs5nHvjApZ99M+dKHrRydd2cPJEL1raWtB97HD8e2bbARX+dzdVGqUnhMjpVSh2BxFVQogOrVeiqAydGcKCajVirvBXYGjgtYz5ikfBqrUPQ1ESN46u557HqrUPF6bgRWhgYCh+dsLvr8DgwFDG/GzTAGBiPIRDB7um7CT7TvWhemHA1mUrJvlaDx6PgtVrVyVsB4HqAPyVfkSjMYRCEwhUs94Be9v+wkXVGB1Rg7rRkSAC1QHTdUPmmmo+hvNvXAQAjP36Chrv/WhC/pIPfBjDP/0JAODsT3+MxR/4EIZ/+hPc97574J3twX3vuyf+fUovX9vB8pXL4genlfP9iMViptsBOfO7m+q3mFITOf4rFDuDiLD2/x4AHUKIKhvnlXexWAxerweAOvwiEo5YytcFR4NQFAX+Sr+9BS5isWgMXp8XAKAoCsLJdW2Sn22afrC0dt0adB7sSugROnb0RSxfucz25SsW+VoPqQRHg2je3Iza+lqblqD42Nn29YOl5s3NUDzcB1nlm+3B1VgYABC+HsU8z50J+XrgAKgBhW+2Bz+6rAYNPU+9gPCNKH7+6ysFLXOxsmPfM3hmELX1ddwOUij0726632JKrWR7IqSUY4ZhTH0AGpM/I4RYL4QYFkIMf+PoN+wqSk4URUEkEgWQuDFlm6/reu7LCIVCaN/1DC4FL6H3xLfsLXgRMgZhsVgMvuS6NsnPNs1f6Y93ndY11MXPkoTGQ/D7Kwq1iEUhX+shlUB1IGGYAdnb9k+e6EX1QrXOfT4vh3FYFL4RxVzFBwDwzfHgrei1hPzeH3wXKz76CRz+wh5U/eF8hG9Ese7jf4Vzb1zA0i99Cm9Fr3E4U5byve9p29mOp3ftgMejcDtIodC/u6l+iym9GRA5vQpXPpskDWG6G8BY8meklEeklIullIvXrFtjV1FyUttQi0uj6rjh0Yuj8aFL2ebrDn/9a9h/aB82Pv5FLAgswLKVf2lvwYtQfX1tfIz2xYuj8a7RdPnZphm7TcevTAYOwdEgFt6z0PZlKyb5Wg9mjOvB5/MiFJqwYxGKjp1tPxaLJUxLP+lB2ek7/33UvE+9DmLJBz5sOjRp6Zc+Fb8m4sQ//wN8sz0J+XoQQunlc99z6GAX1jyyKt7jwO3AXKF/d1P9FlN6QszI6VUotg5nEkL0CCG6AYSllP0Zv+EidfW1iIQjaN/1TPw9ADRv3pY2v+vg8wiFJtB18HnEojGTKVOyuoY6hMOR+N176hrqAEzeUcMsP9s0RVGwdVMz2na2Q9HOgABAKDTBceFJ8rUeAPWHPBRSx79GozEoioK2ne3YuqkZIxdGOIxMY2fbX712FQbPDE6p8+R1Q+ZO/8sPMM9zZ/wOTKf/5QcAJu++5J3tweEv7MHhL+zB8X/+NiI3oth/8jCa7v0YTjz1ApZ84MM4+sr/5Vj5i0m+toMXj72EoYEh7N7Vjg1rH0Xfqb6U20GpK/TvbqrfYkrP7ddECCllwWaWzs3fX3dHQUpQ+R1z8M6tm04Xo2TNLCtn/TuEde+smWXlEB+f73QxSpJ8Rb1eg+3fGdz3OGtmWbm7752q6Rn7P3M6Nl5R9amCLB8fNkdERERE5DKF7FXIBYMIIiIiIiKXKeSdlnLBIIKIiIiIyGXc/sRqBhFERERERC7DnggiIiIiIrJE2HoT1dvHIIKIiIiIyGXYE0FERERERJa4/e5M7u4nISIiIiIi12FPBBERERGRy8zgcCYiIiIiIrLC7cOZGEQQEREREbkML6wmIiIiIiJLeItXIiIiIiKyhD0RRERERERkyQxeE0FERERERFawJyJL5XfMcboIJW1mWbnTRShprH/nsO6dJV+54nQRShrbv3NY95QJ786UpfDv/qfTRShZvnf9Ad65ddPpYpSsmWXlePvWDaeLUZJmlc1m23cQ275zZpXNBgCI1QscLklpki9ews3fX3e6GCWrWE5csyeCiIiIiIgs4d2ZiIiIiIjIEj6xmoiIiIiILOE1EUREREREZAmviSAiIiIiIkvYE0FERERERJawJ4KIiIiIiCyZ4fK7M7m7dERERERE5DrsiSAiIiIichkOZyIiIiIiIkt4YTUREREREVnCnggiIiIiIrKEPRFERERERGQJgwgiIiIiIrKmlIczCSGWA2jS3m6XUobtnF++7Wndh2gkigr/e/DXj38+q/zTr/4TXhs4CwD4682fg+JREIvG8OyuDux5rr2g5S9GbTvbEY1EUFHpx+YtG7PKT5UGAIpHiadFozG0tbThQOf+Ai1N8dm9sx2RSBT+ygpsMql/s/zktOBoEM2bt2NBYAEAoKXtaXg8ilb/u3Ggc19Bl8nt8tnmk9P6TvVh4MwQAGDTlo3weBTTbYMS5WM7AICtm7YBALxeD55u3VG4BZgmDj/cirmzvRj77Tie6Dlomg8A4ZvReL53loKjq3dj5dc2FbKo00b7rmcQiUTg9/uxccsXs8qPRWNo27kb+w+p+/bgaBDbHjf8BrQ+DcWjFG4hphG390TY9pwIIUQVgIeklBu0V9iuednhtYGz8FdWYM9z7fD6vLg0+nrG/InxCfT/42k8uXMbnty5LR5A/P3J/xevB19PMSfSDZ4ZROV8Pw507ofP50VwNJgx3ywtOBrEmkdWoaV1BybGQwiOBhGNxtDb04tg8JJDS+d+g2eG4J/vx4HOffCa1v/U/FTfaVz6AA507sOBzn3xAKK3pxeXWP8J8tXmzdJC4yGceuVVtLTuQEvrDng8ium2QYnytR2cPNGLuoZaHOjch8DCAPpO9Tu0RMXpgep6XP7NFaz82ia8dT2Me++qTshftvhBnPvFCB59eScA4N67quGdpWD9n6xETdJnKTuDA0PwV/qx/1CKtm+SH4vG0NvzrSn79samRuw/tA/7D+1jAHEbhBA5vQrFzofNLQdwVQjRLYTos3E+tnht4Cw+GPgAAOA9/vfg7OBwxvz/3vc9eLwe7Gndhy+sU89EKR4Fn1nzKcxRuBFlMjAwhIB25sLvr8DgwFDGfLO0QHUA/ko/otEYQqEJBKoD8HgUrF67CgrXQ0qDA4MIBAIAAL/fP6X+zfJTfSc0PoHOg13xH6HJ+p9TqMUpCvlq82Zpfa/2w+P1om1nOzasfRQATLcNSpSv7WDhomoER9QTGMGRIALVCwq7IEWuaWEDzr85CgAY+20IjdX1CflL3rsIw7+4CAA4+/OLWPxHixB5O4b9rxxD+Gas4OWdDobODGGB1k4r/BUYGngtY77iUbBq7cNTfltDoRC6Dj7PExW3SeT4r1DsDCLmAbgspdwAoE8I0WjjvPIuFrsOj9cDAFCUOYiEIxnzI+EI/JUVeHLnNtxff398WBNlJxaNwevzAgAURUE4uc5N8lN9Rx1S04za+toCLkFxU+tSb9PK1DZvkm+Wpv+YrFm3Bp0Hn0doPFTApSgu+WrzZmnhcASV8/1oad2BuoY6DJ4ZBMBtI5N8bQd6gNa8We2V9lf6C7gUxc9X7sHVG2rdh29GMW+OLyFfDxwANaDwlfME0e2KxWLw6sc1HpO2nyFfN/kbsBpdz32ZvwG3oZSDiLcAnNf+HgNQk/wBIcR6IcSwEGL4pa//jY1FsU5R5iAaiQJQAwb9BzpdvtfnTeideD3408IWusgZd0qxWAy+5Do3yU/1nUB1AN3HDgNA/OCJ0lPrUm/Tsalt3iTfLM1f6Y8PY6prqJ1yJpcm5avNm6X5fN6E3gl9KB+3jfTytR2cPNGLwMIAuo+9AK/Py+FMFoVvRjF3tlr3vnIP3roeTsjvHT6FFUsexOGHW1H17kr2PuSBoiiI6Mc1UZO2nyFfpw95UjwK6uprp/RoUPbsGs4khFguhGgUQqw3yfMJIWq0z3Skm46dQcR5TAYOVVADiQRSyiNSysVSysWrHvmsjUWx7v76JfEgIHgxGA8O0uV/MPCBeNovQ7/Ee/zvKWyhi1x9fW38QOfixdH4AVC6fLM0Y/epz+dFKDRRoCUobnX1dQgG1bobuTgSH56RLt8szVj/oSsh+P08A5tKvtq8WVogsCCeFgpNqIEEt42M8rUdxGKJB7WRiPlZWzLXN3Imfm3Dkvcuig9tMlp6YG38mogTr71S0PJNR7UNtbg0qu4zRi+OxocuZZuvS/gNGA+horLCphJPf3b0RGg3PYKUsl97nzxSaCWAxVLKk1r+lEBDZ1sQoRXubiFED4C79cIUi/vrlyASjmBP6774ewB48vEdKfPvr1+C0PgEnnx8B0LjE3hg6Z8CAL7y3Nfwy9Av8ZXnvoZYlGdLUqlrqEM4HInfPaauoQ4AsHVTc8p8szRFUe9As3VTM0YujGD5ymUAgEMHuxAKhXDoYBeiXA9T1DXUIhKOYHe8LtXhLvodZszyzdIURcHWTduwe2e7eiZKm07nwS6EQuq1Eqx/Vb7afKq08SshbN3UjPErITQ92JRy26BJ+doOVq9dhcEzQ9i6aRtGLoyyri06PTqAeXN88TswnR4dAACc+HwnAPUuTIcfbsXhh1tx/LVXEHlb3afsXbEFVe/2Y++KLfDO4hAnK+rq1XbcvuuZ+HtAHZKXLr/r4PMIhSbQdfB5xKIxKIqC5s3b0L7rGSiKEv8cWWfTcKYlmDyxP2WkkHaC/4j21rQTIF4+KWWOi5Zf4d/9T3cUpAT53vUHeOfWTaeLUbJmlpXj7Vs3nC5GSZpVNptt30Fs+86ZVTYbACBW84JvJ8gXL+Hm7687XYySVX7HHHffO1Vz8dr5nI6N75l73wYAxh6EeGAghOgG0C2lPK/1QjRJKbcnT0O7y+p27dpmU3zYHBERERGRy+R6kbQWMBxJkR0GMDeLySxPF0AA9l4TQUREREREObDpwuqzAHza31UApjyGQQixXEq5T/t7yo2RdAwiiIiIiIhKgHaNcpU2lMlnuMC6T/u/EUCHEOKcEOIc0vRacDgTEREREZHL2PXMB72XAUC/Ia1J+78fwN3ZTIdBBBERERGRyxTywXG5YBBBREREROQy2Tw4zkkMIoiIiIiIXIY9EUREREREZAmDCCIiIiIisoTDmYiIiIiIyCIGEUREREREZAF7IoiIiIiIyBJeE0FERERERJYwiCAiIiIiIks4nImIiIiIiCxxe0+EkFI6XQadawpCRERERNOWu4/ONVeuX87p2Hj+nLsLsnyu6Yl459ZNp4tQsmaWlbP+HcT6dw7r3lmsf+fMLCsHwN9ep8wsK4do8jtdjJIl+0JOFyErHM5ERERERESWuH04E4MIIiIiIiKXcXtPxAynC0BERERERMWFPRFERERERC7D4UxERERERGQRgwgiIiIiIrLA3SEEgwgiIiIiItdx+4XVDCKIiIiIiFyHQQQREREREVng7hCCQQQRERERkQu5O4xgEEFERERE5DJuvyaCD5sjIiIiIiJL2BNBREREROQyJfuwOSHEegArtLdzARyXUu6za352aNvZjmgkgopKPzZv2ZhVvlna1k3NiMViqPD70dK6o6DLUCzyVddmadFoDG0tbTjQuT9hegCgeBTT+ZUqO9s8twNzdrZ9s3bOtm/uxWMvYeTCSMJ+Qme1zpLXRXA0iK2btyEQWAAAaGlrgcej5HcBily+tgOz/T33Pbk5vHEv5np8GPvVFTzx9WdN8wEgfCMazzdLo9y5PYiwbTiTlPKIlLJJStkEYLjYAojBM4OonO/Hgc798Pm8CI4GM+abpZ080YuWthZ0Hzsc/x4lylddm6VFozH09vQiGLwUn15wNIg1j6xCS+sOTIyHpsyvVNnZ5rkdmLOz7Zu1c7Z9c9FoDONXQqZ5Vuss1TptXNqIA537caBzPwOIJPnaDsz299z35OaBez+Ky796Eyt3P4q3otdw7/sWJeQv+9gncO5nF/Bo1xMAgHvft8g0jaY326+JEEIsB3DO7vnk28DAUPyskd9fgcGBoYz5ZmnLVy6L/2BUzvcjFosVcCmKQ77q2izN41Gweu0qKMrkj3agOgB/pR/RaAyh0AQC1YFCLKbr2dnmuR2Ys7Ptm7Vztn1znQe7sOaRVaZ5qeps66ZmbN3UjJMnehM+n2qdToyHcOhgFwM3E/naDsz299z35Kbpvv+A8z+7AAAY+9UVNNZ8LCF/yQc/jOGf/hgAcPb1H2PxBz5smka3RwiR06tQCnFh9ZNSyiMFmE9exaIxeH1eAICiKAiHIxnzM31n8MwgauvrClD64pKvus40HaPgaBDNm5tRW19rxyIVpUK1eW4Hk+xu+2btnG0/UXA0CMWjwF/pT/sZY529eOwlPPjxpTjQuR+jI4lBgdm60A9q165bg86DXQiNm/d6lCo79j1muO/Jnm+2B1djYQBA+HoE8zx3JuQbg4QlH/wwfHM8pmk0vdkaRAghqgCM2TkPuygeBRFtpxSLxeDTdlbp8tN9p21nO57etYPd2CbyVdeZpmMUqA6weztJIdo8t4NEdrd9s3bOtp+o82AXJsZDaNvZjmDw0pSeBWBqnY1fCWHgzBAOHeyC4lFw6GAXDh3swskTvabrwl/pjw9jqmuom3KmvdTle99jhvsea8I3opir+AAAvjlevBW9lpDf+/3vYMV/+AQOb9yLqn8/H+HrUdM0uj0ix3+FYndPRA2As6kyhRDrhRDDQojhY0e/YXNRrKmvr42Pq7x4cTTebZouP9V3Dmld5enOdJWyfNV1punojMMJfD4vQqGJvC9TMbK7zXM7mMrOtm/Wztn2p+o+dhgHOvdj05aNCAQWYPnKZQn5ZnVWvTCA+oZabN6yMeG1fOWyjOti/EoIfn9FYRauSORz32OG+x7r+s79M2refw8AtVdBH9pktPSJv4pf/3Dif3w7ZRrdDpHjq0Clk1LaN3EhtgEIZzOc6Z1bN+0rSI4OHexCLBpLuCPH1k3N8bs+mOUnp7147CWceuVVeL1qt97ylcvQ9GCTMwuUwsyycrxz66ajZchHXadL6+3pxbIVy7B23RpEIxF84+svIRpRz2KZ3Y2lkNxQ/zq72nwoNOHK7cANdW9X2w+Nh6a0c7M0J7mh/nXJd/XR10GqOtPv2BSNRKbUo9m66DzYBY/X65q7Ys0sKwcA19R/PrcD4/6+t6fXtfse0eTuoGbvI0/BN9uTcKelE08fxsrdj8I724OOR54CAPSd/z56v/8d0zS3kn0hd9/2SHPtX3+b07Hxnf/m3QVZPluDCCvcGESUCjf9kJci1r9zWPfOYv07x21BRKkphiBiOmMQkR982BwRERERkcsU8k5LuWAQQURERETkOgwiiIiIiIjIAneHEAwiiIiIiIhcyN1hBIMIIiIiIiKXcfs1EYV4YjUREREREU0j7IkgIiIiInKZQj59OhcMIoiIiIiIXIdBBBERERERWeDuEIJBBBERERGR67j9wmoGEURERERErsMggoiIiIiILHB3CMEggoiIiIjIhdwdRjCIICIiIiJyGbdfE8GHzRERERERkSVCSul0GXSuKQgRERERTVvuPsWveefWzZyOjWeWlRdk+dwURBQ1IcR6KeURp8tRilj3zmL9O4v17xzWvbNY/85i/ROHM+XPeqcLUMJY985i/TuL9e8c1r2zWP/OYv2XOAYRRERERERkCYMIIiIiIiKyhEFE/nBcoHNY985i/TuL9e8c1r2zWP/OYv2XOF5YTURERERElrAngoiIiIiILGEQQURERERZE0L8G6fLQM5jEHEbhBB3O12GUiVUi4UQW4UQitPlKTVa/S8UQjwqhPA4XZ5SI4R4rxBij9PlKEVa279PCHFICDHL6fKUIiFEjRBipxCi3OmylBKt7d8vhPg7AP+H0+Uh5zGIyJEQQgA4IYT4j4b3VDhfBvCfAfwBgC8IIeY6XJ5S0w7gAQBVADY4XJZStBLAn+j7HyoM7YTFGQCfBNAmpXzb4SKVHCHEfwKwRErZCuCW0+UpFUKIKgC9ABoBXAHwprMlIjdgEJG7DwMYBfBFAJC8Qr1gtDPfp6WULQAOApgNYI6zpSodQogZAL4rpXwewN8C+JjDRSoZWt0DwOsADgP4jIPFKUW/A3AcwK8A1AshGoQQZQ6XqdRUA5gvhNgJ4G+EEO92ukAlIgzgESnlswC+A+BOZ4tDbsAgInc/llJ+BsA8IUQNwN6IAooB+IEQQkgpfwtgPtRAggpASvm/AfxICDETwKcAhIUQHxNCvMvhok17Usr/rQUSn5RSvgzgbSHEFiEEg+gCkFL+K4AfQh3K8S4AqwFsEUJ4HS1YabkFQAHwPQA9AD4jhPgDR0tUAqSUV6WUV7W3nwXwG4DHPaWOQUTu9A3nJLTeCCoMqfqtlFIKIfwAxqWUQafLVUqklDeh7j/+AUAX1OEdH3e0UCVCC+ImhBA7ANwFtd5vOluqkvIvAP4ZwE8AdADwAFjgZIFKzHkAZVB7416FGlQsdrREJcLQE3oOwHKAozBKHYOIHGk/5JBSfgXA3UKI/wygwtlSlaR/C+CbQoh6IcQang0vHCnlTSnlD6SU56AO84g4XaZSIIS4A8B/BPBLAOsA/H8AFjlaqBIipfxfUsr9UsoxKeUbAN4N4H85Xa4Scg7AONTrIqIA5gK4mv4rlA/6cQ/U3rh/5U016A6nC1DMtKh8JgAv1DOx/8PZEpWk1QDuA3AWQIeU8ncOl6dkCCEqANwN4CEA1wH8yNkSlQYp5e+FEA9IKX8PAEKIowBCDherpGht/x4AfwH12rjXHS1QCdHa/zcBrBRCfArAzwCwJ7qw/g0AH4B/DyDqbFHISXxi9W3SbvH3Liklz8IWmHZGdiWAfzSM1aQCEkJUAohJKcNOl6XUCCFmGM4MUoEJIf4IQJT7Hmdovc4ztd4IKiDtergyKeUNp8tCzmIQQURERERElvCaCCIiIiIisoRBBBERERERWcIggoiIiIiILGEQQUREREREljCIICLKkhDimhCiTwjRo/3fcRvTuqz9X6X/nS/ppieE8AkhzmX4fsbP5PJZIiKaPhhEEBFZs0JKuUJK2QQAQohttzMx7aFld6f7jBCi6nbmQURElG8MIoiIcnccQFMB5tNTgHkQERFljUEEEVHuHgLQpw3p6RNCbBNC9ABqD4UQ4pyW7tPSurW0bgBztbSE4UBCiPXaZ85p0+gAUKNP38q0s6FNo0/7rs+QVZU0Tf3zU+ZNRESl5w6nC0BEVGR6hBD63+ellPu0g+lGAH1SyhVCiBoATVLK+7ShSB1acFElpbwPAIQQK5MnrH1vg+EzPillWAjRaBg+ldO0UzFMdz2ADgAbDHkbtLw+IUQjgKvJ8wawPdt5ERHR9MEggojImhVSyrBJelhKuU/7+yEAc/VeCQBVAMJIHJZ01WQaD0EdIgUASDGfXKdtSgsGGjF1WJZxGj0AagDMM5k3ERGVIAYRRET5kXzgvkdKeVJ/czt3cjKRl2lrvRpHAawDMAZDL4SFeftymTcRERU3XhNBRJQfYcPfx2E4INfO9vcBWKG998H8uoXjUHsajN+D4Tu3M20ziwEcl1Kex9ReBeM0NgDoTzFvIiIqQeyJICLKMynlef1ZElpSn3btxArtIuphqGf+zb7XbbjQuh/qNQf9AE4LIYallBtymXYKJwCcE0I0QQ2CjL0pY9oF1YsB9GuBBpLnDeBIlvMiIqJpREgpnS4DEREREREVEQ5nIiIiIiIiSxhEEBERERGRJQwiiIiIiIjIEgYRRERERERkCYMIIiIiIiKyhEEEERERERFZwiCCiIiIiIgsYRBBRERERESW/P8o/BR4cjomEAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Build the plot\n",
    "plt.figure(figsize=(15,5))\n",
    "sns.heatmap(matrix, annot=True, annot_kws={'size':10},\n",
    "            cmap=plt.cm.Greens, linewidths=0.2)\n",
    "\n",
    "# Add labels to the plot\n",
    "class_names = ['1', '2', '3', '4', '5', '6', '7']\n",
    "tick_marks = np.arange(len(class_names))\n",
    "tick_marks2 = tick_marks + 0.5\n",
    "plt.xticks(tick_marks+ 0.5, class_names, rotation=25, fontsize=10)\n",
    "plt.yticks(tick_marks2, class_names, rotation=0, fontsize=10)\n",
    "plt.xlabel('Predicted label', fontsize=12)\n",
    "plt.ylabel('True label', fontsize=12)\n",
    "plt.title('Confusion Matrix for Random Forest Model - Yoshii River\\nTraining data - Kano\\&Yoshii River, 75\\% training and 25\\% testing, Accuracy score is 97.14\\%')\n",
    "plt.savefig('data/yoshii_river/out_img/predicted/confusion-matrix-kano&yoshii-data.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66fc5212-47a9-4f93-a2c4-e12ad422b53a",
   "metadata": {},
   "source": [
    "**Create raster of predicted data - Kano river**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d776c817-3cb1-4d8c-b16d-b7773fad0e18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add one image for projection and shape reference\n",
    "kano_img = rasterio.open(\"data/kano_river/out_img/class/kano_class.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8e95c361-8c43-456b-9266-0069b37a4ae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create new dataframe for predicted class and its index\n",
    "rfc_class = pd.DataFrame()\n",
    "rfc_class['id'] = kano_index\n",
    "rfc_class['class'] = preds_kano"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e8de4d1e-89f3-4ef5-8813-a887d18c19e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rearrange index into shape of reference raster\n",
    "indx = list(range(kano_img.read().reshape(-1).shape[0]))\n",
    "Index = pd.DataFrame()\n",
    "Index['id'] = indx\n",
    "df1 = Index.set_index('id')\n",
    "df2 = rfc_class.set_index('id')\n",
    "df2 = rfc_class.set_index(df2.index.astype('int64')).drop(columns=['id'])\n",
    "mask = df2.index.isin(df1.index)\n",
    "df1['cluster'] = df2.loc[mask, 'class']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c92380fa-4d12-43de-a293-930f4642827d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape the cluster array\n",
    "array = np.array(df1['cluster'])\n",
    "n_array = array.reshape(kano_img.read().shape)\n",
    "class_array = np.where(np.isnan(n_array), 0, n_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "163a3456-8c10-46ee-a76d-e02008949e9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data dir\n",
    "data_dir = \"data/kano_river/out_img/predicted\"\n",
    "\n",
    "# Output raster\n",
    "out_tif = os.path.join(data_dir, \"kano_predicted.tif\")\n",
    "\n",
    "# Copy the metadata\n",
    "out_meta = kano_img.meta.copy()\n",
    "out_meta\n",
    "\n",
    "# Update the metadata\n",
    "out_meta.update({'driver': 'GTiff',\n",
    "                 'dtype': 'float32',\n",
    "                 'nodata': None,\n",
    "                 'width': kano_img.shape[1],\n",
    "                 'height': kano_img.shape[0],\n",
    "                 'crs': kano_img.crs,\n",
    "                 'count':1,\n",
    "                 'transform': kano_img.transform\n",
    "                })\n",
    "\n",
    "with rasterio.open(out_tif, \"w\", **out_meta) as dest:\n",
    "    dest.write(class_array.astype(np.float32))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "749b3411-70c2-4d13-9f35-707d8a9b35cb",
   "metadata": {},
   "source": [
    "**Create raster of predicted data - Yoshii river**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "d9fc2fd6-119f-4e1b-8e85-79c761a0e9a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add one image for projection and shape reference\n",
    "yoshii_img = rasterio.open(\"data/yoshii_river/out_img/class/yoshii_class.tiff\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "fb45e575-5d1b-4509-91b5-1be5dd4ca574",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create new dataframe for predicted class and its index\n",
    "rfc_class = pd.DataFrame()\n",
    "rfc_class['id'] = yoshii_index\n",
    "rfc_class['class'] = preds_yoshii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "8923fe0b-86a2-4f9e-8b10-7f63937bebbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rearrange index into shape of reference raster\n",
    "indx = list(range(yoshii_img.read().reshape(-1).shape[0]))\n",
    "Index = pd.DataFrame()\n",
    "Index['id'] = indx\n",
    "df1 = Index.set_index('id')\n",
    "df2 = rfc_class.set_index('id')\n",
    "df2 = rfc_class.set_index(df2.index.astype('int64')).drop(columns=['id'])\n",
    "mask = df2.index.isin(df1.index)\n",
    "df1['cluster'] = df2.loc[mask, 'class']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "1946f77b-9a2c-4c0f-b037-354cea35776a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape the cluster array\n",
    "array = np.array(df1['cluster'])\n",
    "n_array = array.reshape(yoshii_img.read().shape)\n",
    "class_array = np.where(np.isnan(n_array), 0, n_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "736e0e7e-3e6e-40a3-bc77-1cdf074c6293",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data dir\n",
    "data_dir = \"data/yoshii_river/out_img/predicted\"\n",
    "\n",
    "# Output raster\n",
    "out_tif = os.path.join(data_dir, \"yoshii_predicted.tif\")\n",
    "\n",
    "# Copy the metadata\n",
    "out_meta = yoshii_img.meta.copy()\n",
    "out_meta\n",
    "\n",
    "# Update the metadata\n",
    "out_meta.update({'driver': 'GTiff',\n",
    "                 'dtype': 'float32',\n",
    "                 'nodata': None,\n",
    "                 'width': yoshii_img.shape[1],\n",
    "                 'height': yoshii_img.shape[0],\n",
    "                 'crs': yoshii_img.crs,\n",
    "                 'count':1,\n",
    "                 'transform': yoshii_img.transform\n",
    "                })\n",
    "\n",
    "with rasterio.open(out_tif, \"w\", **out_meta) as dest:\n",
    "    dest.write(class_array.astype(np.float32))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}